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Predicting Response of a Proposed Hydraulic 
Control System Using Bond Graphs 


Dynamic response is an important criterion of quality for many hydraulic control 
syslems. It is suggested that power flow modeling procedures, and in particular the 
recent development of power bond graph techniques, provide the designer-analyst of 
hydraulic control systems with a particularly relevant means for investigating dynamic 
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performance as he designs a proposed system. 
followed by digital simulation of the model is applied to a hydraulic system proposed 
for a particular task. 


The bond graph modeling technique 


Predicted response and subsequently measured response are 


given, compared, and discussed. Only generally available data and parameter assess- 
ment procedures were used for the prediction. 


Introduction 


Modern hydraulic control systems work at high pressures—-up 
to 21 MPa (3000 psi) for industrial and mobile equipment, to 
35 MPa for aircraft systems, and over 35 MPa for some situa- 
tions. The components used are precision finished and expen- 
sive. Hydraulic control systems are commonly used where 
dynamic response is important, be it a manual or automatic 
control situation. These factors require that adequate dynamic 
performance be part of the hydraulic control system design 
process. This is the case for the systems of expensive aircraft, 
How- 
ever, many industrial designers lack the expertise or confidence 


spacecraft, military equipment, and some machine tools. 


to perform this aspect of system design. They rely on past ex- 
perience and after-assembly trial and modification to achieve 
To a significant degree, this 
reluctance to undertake predictive dynamic analyses is due to 


acceptable dynamic performance. 


difficulty in forming a mathematical model of the proposed sys- 
tem, from which the dynamic performance can be predicted with 
confidence. The major difficulty is a fundamental inadequacy 
of the modeling procedures conventionally used for hydraulic 
control systems—transfer functions, block diagrams, structured 
and unstructured sets of equations, etec.—to cope with the in- 
herent nature of hydraulic components coupled to form sys- 


tems. Factors to be considered include: 


1 Hydraulic control systems are usually of small scale, in 
that they consist of relatively few components or subsystems. 


9 


Responses required are often of large scale; typically, a 
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load may be driven from one extreme position to the other. 

3 Nonlinear characteristics are inherent. 

4 Hydraulic systems are made up by coupling more-or-less 
standard components in a building block manner. 

5 Models of components should be reusable. This reduces ef- 
fort and cost of response prediction. 

6 Hydraulic systems are power transfer systems; they work 
via power interchange and transformation along specific paths 
between components. Large resisting loads are usually involved. 

7 The system designer needs a formal procedure whereby 
he can form a mathematical model of his proposed system. 

8 The designer needs to be able to readily manipulate his 
model, simulating changes in the layout of his proposed system, 
in order to satisfy performance specifications. 

9 The model solution procedures available. When only 
manual solutions were possible, coarse linear approximations were 
necessary. Modern digital simulation procedures such as MIMIC 
and CSMP allow highly detailed and nonlinear models to be 
solved to provide time responses. 


The multiport approach to modeling of dynamic systems al- 
lows recognition of these factors. Models are formed on the basis 
of power interchanges between recognizable functional ports on 
the coupled components. In particular the power bond graph 
technique developed by Karnopp and Rosenberg [1]! from earlier 
work by Ezekial and Paynter [2] offers a formal approach to the 
modeling of powered control systems. It allows formation of 
the model structure to proceed, deferring consideration of the 
specific mathematical relationships to be used until finalizing 
the model for solution. The model structure is closely identifiable 
with the functional structure of the proposed hardware system, 


'Numbers in brackets designate References at end of paper. 
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and component models can be connected and disconnected in a 
manner analogous to the physical interconnections of adjacent 
hardware components. The assignment of causality—the desig- 
nation of which variables are to be used to compute other vari- 
a largely separate and formal procedure which pre- 


ables—is 


pares the model for solution. 


The present paper describes the use of bond graphs in a hy- 
draulie control system design situation. The model developed 
and used is modified later to provide information on the effects 
of variations of the proposed model on predicted response. The 
results of actual system response measured subsequent to the re- 


to the predictions. 


sponse predictions made at the design stage are given and compared 


Appendix A provides a table of those bond graph terms and 
symbols most useful in modeling hydraulic control systems. 


The Proposed System 


The requirement was to position a load vertically using a 
hydraulic cylinder to actuate a pivotted beam in a manner similar 
to many mobile equipment applications. The specifications, with 
load and travel compatible with laboratory investigation, were: 





Nomenclature 


an area 

effective area of hydraulic 
actuator (cylinder), head 
end 

effective area of hydraulic 
actuator (cylinder), rod 
end 

effective bulk modulus of oil 

effective bulk modulus of oil 
at pump discharge 

effective bulk modulus of oil 
at the actuator 

effective bulk modulus of oil 
at head end of actuator 

effective bulk modulus of oil 
at rod end of actuator 

a capacitive energy storage 
effect 

hydraulic capacitance of 
pump discharge 

hydraulic capacitance at 
head end of actuator 

hydraulic capacitance at rod 
end of actuator 

general bond graph symbol 
for an effort variable 

a force 

the 

vyl- 


net force exerted by 


actuator (hydraulic 
inder) 

frictional force resisting ac- 
tuator motion 

net force acting at the driven 
load to accelerate it 

force exerted on the driven 
load by the beam 

coulomb friction (stiction) 
contribution to Ffa 

an inertia effect 

effective total inertia of driv- 

(actual load 
plus actuator and beam 
moving parts referred to 
the load) 

a coefficient 

valve orifice flow coefficient 
associated with orifice R; 
and input X7 

valve orifice flow coefficient 
associated with orifice R, 
and input X7 

valve orifice flow coefficient 
associated with orifice R; 
and input X7 


en system 
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r 
Pp 

P, 
(and Pa) 


Py, 
(and Pb) 


Pe 


Qr 
Qu 
Qv2 
Qu; 
Qvs 
QQ: 
(and Qa) 
Q2 


(and Qb) 
AQ, 


AQ: 


valve orifice flow coefficient 
associated with orifice R, 
and input X7 

coefficient of viscous friction 

lever ratio of beam 

a pressure 

pressure at pump discharge 

pressure at valve control 
port A, which is connected 
to head end of actuator 

pressure at valve control 
port B, which is connected 
to rod end of actuator 

pressure at tank port of con- 
trol valve 

pressure drop 
orifice Pz 

pressure drop 
orifice R 

a volumetric flow rate; also 
general bond graph symbol 
for a flow variable 

net pump discharge rate to 
control valve 


across valve 


across valve 


flowrate lost due to compres- 
sion of the oil 

flowrate lost to 
leakage in pump 

jowrate lost through reef 


internal 


valve 
flowrate through valve re- 
sistance R, 
flowrate through 
sistance R; 
flowrate through 
sistance R; 
flowrate through 
sistance R, 
flowrate received at 
end of actuator 


valve 

valve 
ralve re- 

head 


flowrate received at rod end 
of actuator 

flow rate lost due to oil com- 
pression, at actuator head 
end 

flow rate lost due to oil com- 
pression, at actuator rod 
end 

a resistive effect 

hydraulic resistance of pump 
to internal leakage 

hydraulic resistance of pres- 
sure relief valve 

hydraulic resistance of valve 
orifice Ry 


hydraulic resistance of valve 
orifice Rz 

hydraulic resistance of valve 
orifice Fs 

hydraulic resistance of valve 
orifice Ry 

frictional resistance of ac- 
tuator 

a source, i.e., a variable re- 
garded as constant 

constant angular velocity of 
pump 

flow source, related to Sw 
via V0 

gravitational force associated 
with driven load 

denotes “sign of’ 

a power transformation. The 
associated transformation 
coefficient may be shown 
with it. 

a volume 

swept volume of puniwp per 
radian of rotation 

volume of pump dicharge 
side 

volume of actuator nead end 
aide 
lune of actuator rod end 
side 

a displacement 

displacement of actuator rod 

displacement of driven load 

displacement of valve spool 
from its central position 

the port opening of the open- 
center control valve in the 
central position; i.e. the 
underlap; assumed equal 
for all 4 orifices; referred 
to input Xi 

displacement of valve con- 
trol knob from its neutral 
position 

A few other suffices are used and ex- 
plained locally. However they do not 
appear in the final bond graph structure 
or equations, 


Units 

SI units are used. Some descriptive values 
are given also in British units. Mul- 
tiples of base SI units may be used in 
descriptive modes. In equations, only the 
base units are used, 
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Fig.1 Hydraulic load positioning system 


The load to be moved and positioned vertically, 125 kg 
(275 lb). 

Total load travel, 0.6 m (2 ft). 

Time for load to change position following a command 
input, 1 s 

Control to be manual. 

For economy, the smallest pump and cylinder should be 
used, compatible with acceleration occupying about 15 per- 
cent of piston stroke time to limit maximum acceleration. 
Standard hydraulic components to be used. 

Maximum system pressure to be 13.8 MPa (2000 psi). 


Fig. 1 illustrates the application, and the hydraulic control sys- 
tem proposed for it. The details of mechanical design and hy- 
draulic component selection, using conventional static relation- 
ships and data are given by Barnard [3]. Conventional estima- 
tions of natural frequencies were used to approximately satisfy 
rise times aud meat) accelerations. Relevant details of the major 
components «re included laier in the paper. The counterbalance 
valve was needed to oold the load »gainst gravity, following a 
deciswon to use an open-cer ter control valve. This type of valve, 
with a high flow gain coeffi ent in the region of null, allows good 
control of speed without. the need for constant pressure supply. 
The literature on dynamic response contains little information on 
the incorporation of this kind of control valve. The pressure 
relief valve in the upper hydraulic line of Fig. 1, set at 6.9 MPa, 
was included to limit acceleration of «he load during the gravity 
assisted downward motion. 


The Proposed Model 


The purpose of the model is to allow prediction of the time re- 
sponse following command inputs at the control valve. The main 
working stroke of the system is lifting of the load. In this case 
the counterbalance valve and the relief valve in the hydraulic lines 
between control valve and cylinder do not affect dynamic re- 
sponse. The model will be developed for this case. 


The Model Structure. The was considered as com- 
prising the four subsystems illustrated in Fig. 2(a). Models 
proposed for the subsystems were developed in recognition that 
primary interest centered on response of the driven load, rather 
than on dynamic performance of the subsystem itself. This 
allowed considerable simplification of model structure. 


system 


@The Pumping Subsystem comprises the electric motor, the 
pump, the pressure relief valve, and the associated hydraulic 
lines up to the control valve. The reservoir and suction line, 
though components in the subsystem, should not affect its dy- 
namic contribution to load motion. The left section of Fig. 2(b) 
shows the proposed model structure. 
valid include: 


Assumptions felt to be 


¢ The pump speed is constant, allowing designation of the 
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input Drive [ 


BO 4 


a. The Hydraulic System 


yy FA reds iim 


| Xa iR = yXm 
Sg 


mT 


Ai Rta 
b. The Bond Graph 
| __ load _ 


Pump and | Cylinder 
-_ r -_ 


Retief Valve | 


Open - Centre 
Controt Valve | 


ad 


Fig.2 Hydraulic system schematic and its bond graph 


pump’s motor as a source of speed Sw. 

¢ The resistance of the short hydraulic line from pump to con- 
trol valve is negligible. 

* The capacitance of the hydraulic line from pump to control 
valve can be lumped with the pump capacitance as Cp. 


Internal pump leakage and pressure relief valve flows are both 
recognized as resistive effects, and are allowed for by the un- 
The 
development. of flow rate by the pump results from the trans- 
formation of speed Sw to flowrate Q via the pump’s volumetric 


specified hydrauli¢ resistances Rip and Rr, respectively. 


dispiacement rate V@ m3/r. 

& more complicated model involving pump anda motor In- 
erioas, friction, and shaft compliance ean be readily proposed 
(reference [3]). If it was felt. that relief valve dynamics affected 
system response, an expanded ~odel allowing for its inertia, 
friction and compliance is veadiiy proposed (reference [3]). How- 
ever, the present purpose is to study resnouse of the load to 
The pumping system 
A .najor advantage of 


input commands at the control valve. 
model proposed was felt to be adequate. ; 
bond graphs is that a component model later found to be in- 
adequate can be expanded (or simplified) with zero or minimal 
effect. on the structure of the rest of the model. 

The location of the half arrows indicating direciions of power 
flow is natural. They can be assigned as the model is developed, 
or after the whole system model is structured and coupled. The 
short vertical bars indicating causality are readily added using 
normal. preference rules (Appendix A), completing the appro- 
priate portion of Fig. 2(b). 

@The Open-Center Control Valve included on Fig. 2 can be 
represented functionally by Figs. 3(@) and 3(b), in which supply 
flow is directed between tank, discharge port A, and discharge 
port B, according to the modulations of the four control orifices 
shown as variable restrictors R;, Re, R;, and Ry The orifices 
open and close in pairs as the spool is moved by command, R; 
and FR, control the flow in discharge port B, and R: and R, the 
flow in port A. A short discharge line to tank justifies neglect of 
power lost to this line from the valve when forming the dynamic 
model. The bond graph model of Fig. 3(c) was proposed for the 
valve, the implication being that the R functions are modulated 
by the input command displacement Xi applied to the valve 
handle. Each R may be a nonlinear function, though it is not 
necessary to decide this at this stage. It is assumed that flow 
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Ga 
e Bond Graph in Circular Form 
(as used in Fia7) 


Fig.3 Development of bond graph structure for open-center control 
valve (each R is modulated by input Xi) 


RL 


forces, hydraulic lock, and other effects on the valve spool are 
minimized by design or dominated by the operator so that they 
do not affect the dynamic performance of the driven load. 
Note that Fig. 3(c) allows for supply flow to go to either pair 
of valve orifices and consequently to either control port. When 
assigning power flow directions, it should be appreciated that 
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(a) Schematic. 


P2\Q2 
Cyl.end1 Ria CyLend 2 


(b) Bond Graph,detailed model 
Fig.4 Detailed cylinder bond graph 


power flows from one of the control ports and into the other. 
Figs. 3(d) and (e) show the bond graph arranged into parallel 
and circular forms which may be preferred. Power flow and 
causality assignment are again straight forwards, yielding the 
appropriate section of Fig. 2(b). The Pumping and Control 
Valve models are readily connected by their common power bond 
yielding the left half of Fig. 2(b). 

@The Actuator (Cylinder) bond graph structure of Fig. 2(b) 
was simplified from the detailed model of Fig. 4 which allows for 
cylinder entry loss (Ra), cylinder end capacitances Ca;, Ca», 
leakage across piston (Ria), piston and rod mass (Ja), piston and 
seal friction (Rfa). Hydraulic power (pressure « flow rate) and 
mechanical power (force « velocity) are related through the ef- 
fective piston areas A; and A», represented by the 7'F elements. 
Ia was made part of load mass Im. Ra and Rla were neglected. 
Note that cylinder capacitance Ca is recognized but not yet 
functionally specified. It can later be specified to allow for hy- 
draulic lines, and for large displacement of the piston. It is 
assumed that severe cavitation will not be encountered, although 
this also can be included in the definition of Ca (reference [3)). 
It is apparent on Fig. 2(b) that the resistances of the short hy- 
draulic lines connecting valve to cylinder have been neglected, 
leaving P; = Pa; Q: = Qa; P2 = Pb; Q. = Qb. 

@The load model includes allowance for the inertias of cyl- 
inder moving parts, beam, and the driven load. Again simplified 
from a detailed bond graph (reference [3]) the proposed model 
(right-hand section of Fig. 2(b)) recognizes the input power as 
force « velocity in its left-hand bond. This power is transformed 
by the beam lever ratio LR to another product of force and 
velocity at the mass. Inertia Jm recognizes the miass of the driven 
load plus the equivalent masses at the load of the beam inertia 
and the cylinder moving parts. The source Sg allows for the 
constant gravitational force resisting upward motion of the 
load, It is apparent on Fig. 2(b) that structural compliance has 
been neglected. 

@The system bond graph structure, Fig. 2(b), is formed by 
direct connection of the four subsystem bond graph structures 
at the obvious common bonds (power ports). 


The Equations. The equation set prepared from the bond 
graph structure of Fig. 2(b) is as follows, the bond graph symbol 
associated with each equation being indicated at the left. 


@ For the Pumping Subsystem 
Cp : Pp = Pp(o) + 1/Cp + [AQp at 


Rip : Qlp 
Rr : Qr 0 


Pp/Rlp 
for Pp < 13.8 MPa 
(Pp — 13.8)/Rr, for Pp > 13.8 MPa 
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0 Junction : AQp = VO+ Sw — Qr — Qlp — Qp 


where AQp is flowrate lost to compression. 

@ For the Control Valve 

R, ° Qu; a 
R, : Qu. 
R; : Qu; 
R, ° Quy; 

0 Junctions 

ieft : Qp 


K, « (X(o) — XZ) « Pl? 
Ke « (X(o) + XZ) « APv,!? ¢ sgn AP, 
K; « (X(o) — X71) « APv;!? « sgn APv, 
Ky + (X(0) + Xi) « P,l? 


Qv2 + Qu; 
Qu, — Qos 
Qn — Qu, 


upper > Qr 
lower : QQ 
1 Junctions 
pper : APoy = Pr — Ps 


lower : APv. = Pp — Pi 


@ For the Cylinder 
Ca, ° P, : P,(o) + 1/Ca; + AQ, dt 
Px(0) + 1/Caz * AQ: dt 


sgn Xa * Fu + Kf + Xa 


Caz ° P, 
Rfa : Ffa 


0 Junctions 


upper A, + Xa — Q: 
Qa - Ai ° Xa 


P, « Ay; — P2* Ar — Ffa 


: AQ: 
: AQ, 


1 Junction : Fa = 


lower 


where AQ denotes flowrate lost to compression 


@} or the Beam and Load 
TF : Xa = Xm/LR 


Fm, = Fa/LR 
Im : Xm = Xm(o) + 1/Im« SFm dt 


1 Junction : Fm = Fm — Sg 


@ Auxiliary Equations 


Input : Xi = to be specified 


Response : Xm = Xm(o) + [Xmdt 


@ Comments 


Most of the equations have well established forms 
Control valve orifice flow is assumed turbulent 

K, to K, are orifice flow coefficients, referred to input Xi 
X(o) is valve underlap, assumed equal for all 4 orifices, 
referred to Xz 

APv, and APv; are the pressure drops across R; and R;, 
respectively 

Tank-side pressure is assumed as zero, leaving P; and P: 
as the respective pressure drops across R; and R2. 

The friction equation associated with Rfa is designated to 
have a Coulomb component (sgn Xa « Fu) and a viscous 
component (Kf * Xa) 

P, and Q; relate to discharge from control valve port A 
P; and Q: relate to discharge from control valve port B 


The Coefficients. Completion of the model requires specifica- 
tion of numerical values for the various coefficients of the equa- 
tion set. 


@ Those following directly from choice of hardware are: 


Sw = 151 r/s (pump speed 1440 rpm) 


ve 1.56 X 10-6 m3/r (0.59 in3/rev) 
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1.14 X 10° m? (1.767 in?) 


0.63 X 10-3 m? (0.982 in?) 


5.1 (average considering angularity; over- 
all length of beam from pivot to load 
center is 1.52 m) 


Im 141 kg 
S¢ 1.38 x 10° WN 


X(o) = 0.019 m 


@ Those readily calculated from manufacturer’s or geometric 
data were 

Cp 2.91 X 10-2 m5/Pa (Volume Vp 2.6 X 10-%/ 
bulk modulus @p 89° x 


106Pa) 


0.505 K 10% Pa/(m*/s) (manufacturer's data} 


© Pa/(m*/s) for Pp < 13.8 MPa (man 


ufacturer’s data) 


0.41 X 10° Pa/(m*/s) for Pp > 13.8 MPa 


@ Those requiring more consideration were 
K, to Ky. The control ports were square. The valve manu- 
facturer provided AP — Q data for each orifice 
for the full open state. Fitting these data to the 
turbulent flow equation yielded 


13.7 X 10-6 m?/(Pa!? « s) 
10.5 X 10-6 m?/(Pa!? « s) 


More detailed AP-Q-Xi data would have been 
appreciated as a system design aid. 


Ca, and Caz. These were taken initially as the ratio of the 
appropriate mean volume and bulk modulus, 
yielding 


= 0.33 X 10-3/535 x 10° = 0.616 K 10-!? m3/Pa 


= 0.251 < 10-3/380 x 106 = 0.664 X 107! m/Pa 
A more complex assessment of Ca will be 
presented in a following section. 


= 120 N (27 lbf; assessed from pub- 
lished data for hydraulic 
cylinder friction, involving 
seal squeeze and pressure 
distribution) 


Kf = 1.57 X 108 Nes/m_ (mean of several reported 
values for hydraulic actu- 


ators). 


The initial values of Pp, Pi, and P: were obtained by giving X7 
an initial value Xi(o) sufficient to generate P;(o) = 6.15 MPa, 
the counterbalance valve setting needed to support the load in a 
stationary mode. The initial value of Xm was taken as zero. 


The Response 


The Input Command. The situation requires load position 
control with position changes to be completed in a specified time. 
The necessary input command Xi(t) at the valve lever is of a 
half wave form. Tests showed that valve stroking motion in this 
kind of control situation approximates a truncated triangular 
wave form with a rise time of 0.06 seconds, and a duration of 0.85 
seconds. This input (Fig. 5) was adopted for the model simula- 
tions. To produce the necessary spool motion, the maximum 
travel of the control knob was 1.9 cm. 
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Valve Displacement 














Piston Velocity 














0.8 


0.6 


TIME, SECOND 


Fig +) Predicted and measured responses. (Full lines show meas- 
ured. 


The Simulation Procedure. The digital simulation procedure 
MIMIC was used on a CDC Cyber 6400 computer. The model 
is described in MIMIC as a set of simple statements describing 
the equations derived from the bond graph structure, and the 
coefficients values assigned. Initial values, modulating input, and 
response variables required as outputs are stated to complete the 
package. 
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Initial Results. Fig. 5 shows the predicted time histories of 
some of the system variables when the desired displacement of 
the load was 0.78 m in 1 second. The corresponding power cyl- 
inder piston displacement is 0.15 m. The full lines on Fig. 5 show 
the actual values of these variables measured when the system 
was assembled subsequent to the simulation. In the experiments, 
pressures were sensed with Ether Type BP6 strain gauge trans- 
ducers. Actuator rod displacement was sensed with a linear 
potentiometer (Penny and Giles Type PGS) mounted in parallel 
with the rod. Velocity of the driven load was sensed with an 
attached linear accelerometer (Kistler servo type) whose output 
was integrated through an analog computer integrator. Spool 
valve motion was sensed using a Schaevitz LVDT attached to 
one end of the spool. The outputs of the sensors were displayed 
and recorded from a 4-beam Techtronics storage CRO, and a 
Sanie UV recorder. Each sensor-recorder channel was calibrated 
in situ. 

The predicted and measured responses agree well in form and 
in time based events. Discrepancies are attributable to approxi- 
mations made in the model form, in the relationships used, and 
the coefficient values chosen. Some of these are investigated later 
in the paper. Also, some aspects of component behavior were not 
foreseen. For example, swashplate angle of the servo-controlled 
pump was reduced transiently on application of the load. This 
resulted in less initial flowrate after the input command than was 
anticipated, and contributes to both the time lag and the dis- 
crepancy in final piston displacement. The two micron filter in 
the pump delivery line caused a noticeable pressure drop; had 
this resistive effect been included in the model, it is probable 
that the measured pressure trace would be closer to the predicted 
characteristic. The capacitance effect of the filter had been in- 
cluded in the evaluation of Cp. It is aslo apparent on Fig. 5 that 
measured pressure P, was slightly affected by premature action 
of the relief valve. Keeping in mind that all parameter values 
used were derived from available component manufacturer’s 
data and from common literature recommendations (typically 
Merritt (4]), the agreement is regarded as quite adequate. It was 
a simple matter to subsequently adjust parameters to bring the 
model and actual responses still closer. 


Effects of Model Variations. The designer needs the simplest 
model adequate to his needs. Experimentation with the model 
can increase the confidence level in his predictions. Bond graph 
modeling allows ready expansion (or contraction) of the model to 
embrace effects or relationships previously neglected. It was 
decided to include: 


1 Allowance for change in hydraulic cylinder volume with 
piston motion, and for other effects mentioned earlier, leading to 
the expressions (reference [5]): 


Ca, = 0.452 X 10-143 + A; « Xa/Ba,; m3/Pa 
Ca, = 0.765 X 10-13 — A; « Xa/Baz m3/Pa 


2 An allowance for beam lever ratio variation during stroke. 
Examination of the lever ratio mechanics led to the approxima- 
tion 


LR = 1,52/(0.288 — 1.81Xa(Xa — 0.203)] 


3 An allowance for load moment arm variation during stroke, 
approximated by 


Lm = L(l — 0.56?) 
where L is distance from load cg to beam pivot center (1.52 m) 


6 is angle of rotation of beam about pivot, measured from 
mid-stroke position. @ is a function of Xa. 


As far as the power bond graph model structure is concerned 


these variations require: no modification for R;—R,: addition of 
Xa modulation to Ca;, Caz, and the 7'F for LR; and the inclusion 
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Fig. 6 Response from modified model. 
simulation results.) 


0.8 SeconD 


(Broken lines are Fig. 5 


of the Xa modulated 7'F Lm in the Sg bond. The model remains 
of 4th order, but is more nonlinear than previously. Fig. 6 com- 
pares some of the results of simulation of the new model with the 
original model results. The differences in predicted performance 
are of little consequence to the system designer. 

The model was further expanded to include separate rep- 
resentation of hydraulic piston mass, beam compliance, hy- 
draulic hose compliance and resistance, and fluid inertia. This 
led to the model of Fig. 7, which is readily developed by expand- 
ing the original bond graph model Fig. 2(b). When the relation- 
ships are compiled, or from the bond graph itself, it is seen to 
be of 10th order with a considerable number of non-linearities. 
The inclusion of the additional effects had little effect on the 
quality of the response prediction for the present system, show- 
ing up mainly as small-amplitude high frequency effects on pres- 
sure traces, 


Conclusion 


The paper demonstrates the application of bond graphs to the 
design of hydraulic control systems. The idea is to ensure at the 
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| Cylinder p Seam and Load-—~ 


Se 2 | 


Expanded bond graph model includes: 
« Capacitance Ch, resistance Rh, inertance th associated with 
line between vaive and cylinder 
« Inertia la of cylin’~- meving parts, separated from im 
« Compliance Cd and “an 3ing Rd associated with beam 
« Allowance for loa®’ mwneat arm change dvring motion, as 
Lm 


Fig. 7 


design stage that a system’s dynamic response will be satisfactory. 
The modular, reusable, locally adjustable, characteristics of bond 
graphs are especially convenieiit in the system design situation, 
where a variety of components and component layouts have to 
be considered in a search for satisfactory system response. The 
clearly defined formal and separate steps in model formation 
using bond graphs are particularly suited to the hydraulic control 
system designer’s needs. 
These steps can be summarized: 


selection of bond graph structures for the various hardware 
components or subsystems proposed for the system. Ef- 
fects which are judged to affect the response of the particular 
system are recognized (as S, R, C, I, TF, effects) without 
much regard to specific relationships or coefficients, These 
will be decided separately and later; 

bonding of the component bond graphs to form the system 
bond graph. This is directly and physically relatable to 
the interconnection of the hardware component’s power 
ports to form the hardware system. 

assigning power flow directions to provide a sign convention 
when it comes to preparing model equations, 

assigning of causality, using established preferences and 
physical sense. This provides an important check that the 
proposed structure is reasonable. Causal conflicts must be 
resolved. 

Preparation of the equation set. The bond graph structure, 
complete with power flow direction and causality, provides 
the skeleton from which a complete and compatible set of 
equations can be prepared in a form suitable for simulation. 
Precise equation forms are decided at this stage, although 
they can be modified later if required. 

Decision of the numerical values of the coefficients (and 
possibly functions) in the equation set. 


The work was undertaken in simulation of the hydraulic sys- 
tem design situation, in that only common components, equations 
of accepted forms, and numerical values available from manu- 
facturer’s data or in common references were utilized. The bond 
graph structure provided the key to what numerical data were 
required, and in what form. The component bond graphs utilized 
are simplifications of what would be required to study in detail 
the dynamics of each of the components. In the present case, as 
in the case of most relevant design situations, the primary in- 
terest lay in the response (motion) of the driven load, and in the 
pressure developed in the hydraulic cylinder. 

The initial prediction of system response correlated well with 
the response measured from the subsequently constructed hard- 
ware system. The bond graph structure was later adjusted to 
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investigate the sensitivity of response to various model refine- 
ments. These confirmed the validity of the fairly wide range of 
simplifying assumptions made in forming the original bond graph 
structure and its resultant equations. 

It is considered that bond graph modeling coupled with digital 
simulation provides the hydraulic control system designer with a 
superior technique for predicting, and improving if necessary, 
the dynamic response of a proposed system. There is a quite 
striking affinity between the modular and reusable natures of 
bond graphs and of hydraulic control components. 
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APPENDIX A 
Common Bond Graph Terms and Symbols 


TERM SYMBOL | AS USED IN PRG 


MEANING AND COMMENT 





Power Bond 


Effort 


Flow 


Power flow in the hond equals the product of E & Q 
Effort means potential 


Flow means flowrate 





Source 


Effort is constant 


Flow is constant 





Transformer 


El transforms to E2, Ol transforms to Q2 such 
that E1.Ql = E2.Q2 


El transforms to Q2, 91 transforms to E2, such 
that El1.91 = E2.02 





Resistive Effect 
Capacitive Effect 


Inertive Effect 


Recognises a resistive effect 


Recognises a capacitive effect 


Recognises an Inertia 








O-Junction 





l-Junction 


Algebraic summing of flow variables <t same E 


ice. QL + Q2 + Q3= 9 


Algebraic summing of effort variables at same 0 


i.e. El + E2 + E3= 0 





Power Flow Direction| Gam 
(half 
arrow on 
one end 
lof each 


bond) 


Shows direction of power flow in each bond 


Easily decided for main power flow sequence 


Power flows into R, C, I, elements (a convention) 


Gives sign convention for equations to be formed 
for junctions 





Causality 
(short 
bar ac- 
cross 
one end 
of each 
power 
bond) 
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Shows which of che 2 variables on a bond is 

cause and wich is effect 

Decides how R, C, I, equations are to be arranged 
Only 1 causal bar can be at an O-junction 


Only 1 bond can be at a l-junction without a 
causal bar 


C elements have causal bar away from the C end 


ii bonds have causal bar at the I end 





Effort is directed at the causal bar, 
flow away from it 
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A need exists for inexpensive interactive simulation of nonlinear bond graphs. 
shown how existing block diagram oriented simulation languages can easily be adapted 
to accept bond graphs freely mixed with block diagrams. 
writing the input structure table and for remedying some often occurring causal con- 
flicts and algebraic loops. 


Simulation of Bond Graphs on Minicomputers 


It is 
Procedures are stated for 


Examples are given of hands-on simulation using the 


THTSIM language, written for PDP-11 minicomputers. 


Introduction 


Bond graphs are a clear and systematic tool for performing the 
difficult step of setting up a model of dynamical physical systems 
[1, 2].!. When the bond graph model has been obtained important 
properties of the system and of its computability can already be 
recognized. 

In a next step system state equations, ordinary differential 
equations or transfer functions can be obtained in an algorithmic 
way. In education and research, however, it is also a general desire 
to be able to study the time response of a system by simulation. 
For linear systems the ENPORT-program [3] fulfills this need 
in an elegant way, because it accepts the bond graph as input 
data. It presents the analyst, apart from the time response, with 
the system matrix and its eigenvalues. The linear approach has, 
moreover, the advantage that algebraic parts of the system are 
easily handled. The required linearity is on the other hand one 
of the restrictions to the general use of ENPORT, another re- 
striction being the requirement of a large digital computer and 
the not too easy implementation of the program. 

To bypass these restrictions the analyst may use a general 
simulation language like CSSL, CSMP or MIMIC, which most 
computing centers have at the disposal of their clients. In that 
case it will first be necessary to derive the system state equations 
from the bond ,, aph. 

Not too much programming experience is needed to write 
these equations in the simulation language format. This is 
generally FORTRAN, extended with many easy-to-use dynamic 
functions and output possibilities. There is one common dis- 
advantage to the use of ENPORT and of the general simulation 
language: the analyst never has an immediate contact with the 
response of the system. A structure or parameter change in the 
model may even take up to several hours of turn-around time. 


‘Numbers in brackets designate References at end of paper. 


Contributed by the Automatic Control Division for publication in the 
JOURNAL or Dynamic Systems, MEASUREMENT, AND ConTROL. Manuscript 
received at ASME Headquarters, December 8, 1976. 
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Hands-on Simulation on Minicomputers 


It has for long been recognized that a minicomputer can be 
used for digital simulation without these disadvantages. The 
user has an immediate contact with the simulation result, which 
can be displayed on a screen. He can stop the simulation, change 
parameters and restart at will. When experimenting is done the 
model can be punched in papertape or stored on disk or cassette. 
On returning to the machine the stored model can be reexecuted 
within seconds. The main disadvantage to this kind of con- 
versational simulation has been investment cost. The large de- 
crease in price of minicomputers however has outdated this 
argument and microcomputer technology will even more con- 
tinue to do so. In spring 1977 investment cost was between 
$12,000 and $25,000, the highest price including a mass memory 
and a graphic display with hard-copy unit. 

Many simulation languages for use on minicomputers are 
available [4, 5, 6] one of them being THTSIM, developed by 
Kraan and Meerman at Twente University of Technology 
(THT) [7] for the PDP-11 series. The main characteristics of 
THTSIM? are: 


¢ Model specification in block diagram form 

¢ Floating point arithmetic (no scaling required) 
Over 40 algebraic, logic and dynamic function blocks (in- 
cluding non-linear functions, PID-controller, sample-hold, 
noise generator) 
Each integrator can be Euler, Adams Bashforth II or out- 
put-limited 
Simulation can be executed in real time (rather slow sys- 
tems only, smallest time constant > 0.1 sec) 
Analogue input and output to external apparatus possible 
The program is written in assembler and occupies 8k 16-bit 
words 
The model requires 4k 16-bit memory space for every 400 
function blocks used 


*THTSIM is under certain restrictions available via the author, THT, 
P. O. Box 217, Enschede, The Netherlands. 
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The essential characteristics of THTSIM and of many other 
minicomputer simulation languages is the block diagram input 
format. This paper describes how these languages can be simply 
adapted to accept input data in the form of bond graphs, freely 
mixed with function blocks. 

Because a causal bond graph can be considered to represent a 
block diagram the input specification for the latter shall first be 
considered in some detail. 


Block Diagram Input Specification 


From the definitions of some typical function blocks (Fig. 1) 
can be seen, that every block has one output variable (y). Fur- 
thermore blocks may have several input variables (x). Param- 
eters, like p or y(0), are indicated under or above the symbol. 


function 





summation 


1 





multiplication by a constant 


ypix 





multiplication of variables 


yellsx 





integration, Adams Bashforth |! 





T : integration time step 


Definitions of some THTSIM function biocks 


Fig.2 Block diagram of RC-network 


output inputs 
1 


2 


GAl 
INT 
GAI 


GAI 


Fig. 3 Structure table input 
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Fig. 2 shows a model in block diagram form of an RC-network 
in which some of the above defined function blocks appear. This 
model is prepared for simulation by assigning a number to the 
output variable of every block. As each block has only one output 
variable this can also be interpreted as each block having an 
individual number. The structure of the above model can be 
completely stated in the form of a structure table (Fig. 3). In 
this form programs like THTSIM accept the structure of the 
model. 

It can be remarked that each line essentially represents an 
elementary equation of the model: 


output = function (input) 
For example, for the integrator 5 
ys = INT (ys — yz) 


After the structure input has been typed in, the program asks 
for other necessary data, like parameters, plot variables and 
timing data. After that, commands can be given to display or to 
record a plot, to print an output table, to change the structure or 
the parameters, to store the model on disk, etc. 


Adapting the Simulation Program to Accept Bond 
Graphs 


A causal bond graph is a compact representation of a block 
diagram. The causal stroke at the bond(s) of each element indi- 
cates whether in the computational procedure the e-variable is 
the input and the f-variable the output or whether the inverse 
is true. So for each element the input and the output variable 
are known. The bond graph element itself behaves in the sim- 
ulation procedure like a functional block, that performs an opera- 
tion on its input resulting in an output. Consider for example an 
R-element, Fig. 4. 


" f [oa | 
‘ie 9g 


Fig.4 R-element and associated function block 


With the indicated causality e is the output variable, f the 
input and the function performed is the multiplication with a 
constant. In THTSIM like in many other simulation programs 
this function would be called GAIN. Adapting the program to 
bond graph input can in this case be done by introducing a seem- 
ingly new block, named R, which performs the same function 
as GAIN. In the structure table, the R-element will then appear 
in the form 


output = R (input) 


As Table 1 shows, most bond graph elements can be imple- 
mented in the program by using new names for existing function 
blocks. Adapting the program to accept a storage element how- 
ever requires the definition of 1 new composite function block 
named C or J, Fig. 5. This block is a combination of integration, 





Fig.5 C-element and associated composite function block 
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division by a constant and addition of an initial condition. When 
used this block appears in the structure table like 


output = C (input) 


In the parameter list, two parameters are needed, i.e., the values 
of C and of e(0). It will be clear that an J-element requires the 
same composite function block, now bearing the name J. 
Storage elements with derivative causality require the definition 
of another new composite function block named DC or DI as 


shown in Table 1. The D stands for “Derivative” or “De- 


Table 1 Definition of bond graph blocks 





bond graph bond graph 
element | block 


equivalent 
function (block (s)) 


| 
} 














§ 





—} + 


oO 





sign(De). (Del G) 





; 


fe 
+ 





t 
7 


if 


f 


T: integration time step 
Initial cond’s De(o) and f(o) 





’ + 2 ist — vt 
a= bay + Ct, —2t9 


oc 
“={o1 } 


T: integration time step 
initial cond’ s: Df(o) and e(o) 





t 
7 
nn 
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pendent.”’ The differentiation procedure is the inverse of the 
2nd-order Adams-Bashforth integration, made realizable by an 
extra unit time delay. Certain precautions apply to the use of 
these elements, as explained in a later section. 

R-elements can have two causalities, both of which are realized 
by a GAIN-block. 
the author has found it useful to distinguish between both cases. 


To avoid errors in assigning the parameters, 


This is done by defining a G-element apart from the R-element. 
Mechanical and hydraulic R-elements often have a nonlinear 
“absolute square” For this purpose RSQ and 
GSQ have been defined as shown in Table 1. T7'ransformers and 
gyrators have two ports with a distinct output variable. Two 
function blocks are needed as shown in the next section. 

To increase the readability of the structure table, SE and SF 
blocks can be used. They name the links between the bond graph 
and other parts of the model. 

Junctions constrain both efforts and flows in such away, that 
they specify one kind of relations between e-variables and a dif- 
ferent kind between the f-variables. So, in a block diagram two 
kinds of operations are involved, Fig. 6. 


characteristic. 


fi+fo=f; 
f; 


fe 
C2 


Fig.6 0-junction and corresponding operations 


For a zero junction the summing of the f-variables is from a 
computational point of view the important operation. The 
identity of all e-variables is in a given structure more a definition 
than an operation, but it is on the other hand very essentially 
linked with the concept of the 0-junction. 

A practical solution to this ambiguity can be found. Junctions 
are eliminated from the structure by considering them a part of 
their causal-dominant element. The result is, that the summing 
operation of the junction is becoming an integral part of the 
operation of the element. The defining operation of the junction 
will spread the causal output of the element to all bonds con- 
nected to the element-junction pair. 


f, e @3= R.C f:+f2 
3 
pik =RSt 


— fs 


Fig.7 Typical bond graph element with dependent junction and their 
block realization 


In Fig. 7 an R-element dominates the e-variable of the 0- 
junction. This implies that all variables of the junction are 
named ¢;, after the output variable of R. Furthermore the 
summing operating =f has become an integral part of the op- 
eration of the function block. In Table 1 all bond graph blocks 
are defined in this way as indicated by the Ye or Lf of the input 
variables. 

In some cases junctions cannot be eliminated, because they 
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are causal dependent on more than one bond graph element, via 
another complementary junction. For this case a 0-block and 
a 1-block are defined, as shown in Table 1. Examples of such in- 
dependent junction pairs occur in the next sections. 

As is shown in the examples the elimination of element de- 
pendent junctions as separate elements results in a very clear 
structure table. Its readability is comparable to a well ordered 
set of equations. This makes the structure table also very useful 
for the case that an expression-oriented simulation language like 
CSMP is used to simulate the bond graph. 


Preparing the Structure Table From a Bond Graph 


The structure table, which specifies the bond graph to the 
simulation program, can be obtained in a systematic way: 


1 The bond graph is augmented with sign convention half 
arrows and causal strokes. 

2 Causal conflicts and causal uncertainties are removed as 
discussed in a later section. 

3 It is indicated to which elements the junctions belong. 


Explanation: Nearly every junction is a summing point for 
the input of one specific element. This dependency can be in- 
dicated by a dotted line around the element-junction pair. Ex- 
amples (Fig. 8): 


Cc I G: 
i | | i. 
bOI ah 0 TE 


Fig.8 Element dependent junctions 


0- or 1-junctions which are not a part of an element will occur in 
pairs, called an independent junction pair. (This corresponds to a 
“junction causal bond”’ in reference [8].) They are indicated by a 
dotted line around the pair, Fig. 9. 


io} I — 


Fig.9 Independent junction pair 


4 Output variables of all bond graph elements are indicated 
by a number. For readability of the structure table this number 
can be preceded by the letter e or f of the variable type. Inde- 
pendent junction pairs, transformers and gyrators have two dif- 
ferent outputs to be indicated by two numbers. 

5 The structure table is written and can be typed in to the 
computer. The structure table is obtained as follows: For each 
output variable the associated bond graph element is named, 
followed by its input variables. The latter have a positive sign 
if the half arrow points towards the element. 

6 Prior to simulation a check of the initial conditions and a 
check of all loop gains is advisable: Signal loop gains are directly 
written from the bond graph, reference [9]. They are checked to 
see whether unwanted positive feed back loops exist, generally 
indicating an orientation error. From the highest loop gain the 
maximum integration step size is estimated. 
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Example 1 RC-Network 


ideal physical models 





bond graph 


fully augmented bond graph. 


f 
G,1 pair. 


is output of the encircled 


Inputs to the pair are both 
bonds connected to the pair, i.e. 
e) and ~€3- 

STRUCTURE TABLE THTSIM structure input 
output, elem. inputs 


» & a % “2 





e 
f 


S « e, ~e3 
S f, -f 
G 


1 

2 

_ * Q 

4 3 


» &3 


Fig. 10 Deriving the bond graph simulation structure table for an 
RC-network as in Figs. 2 and 3. 


Example 2 Propeller Drive System 


Ri 


Fig.11 Propeller drive system incorporating a torque controlled DC- 
motor, vee-belt transmission (C,R2) and nonlinear (absolute square) 
load (R3) 


Fig. 12 Combined bond graph-biock diagram. An extra 0-junction 
is introduced at the propeller for extracting the load torque 


‘ —— 4 


Fig. 13 Fully augmented bond graph. e8-/9 is an independent junc- 
tion pair 
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Structure Table 


Output Element Inputs 
é1 SE 12 

fe ’ G a~— & 
€3 GY 

4 GY 

fs I 

& © 

er R 

3 9 01 
fo » a8 
fio ’ I 

eu 9 RSQ 
12 ‘ PID 


13 CON ; (desired value) 


NOTE: 'The 0 and 1 elements may be interchanged, as both elements perform 
@ summing operation. 


PROBLEM: EXAMPLE 2 PROPELLOR DRIVE 
L 

TIMING? 

RANGE: 

PL8SLKS: 


parameters - - block 





THTSEN 


FULL-LOAD 


OU-LOAD 





—_ 


ed 


(b) 


Fig. 14 THTSIM listing and X Y-recorder output 








The procedure for the modeling and simulation of a nonlinear 
propeller drive is shown in Figs. 11-14. The computatien time 
for the shown responses on a setpoint change was 25 seconds on 
a PDP 11/10 without hardware floating point arithmetic. 
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The Removal of Causal Conflicts and Causal Un- 
certainties in Some Often Occurring Cases 


In the foregoing we have restricted ourselves to the simulation 
of bondgraphs using conventional simulation programs. This 
limitation excludes systems with computational problems that 
require special techniques. For example, systems incorporating 
a nonlinear algebraic part (i.e., a resistance network) can not be 
solved; if one is willing to solve a stiff system, though, one may 
introduce small capacitances in the network. 

Bearing in mind the foregoing restriction, it will be shown by 
examples that some often-occurring causal problems can be re- 
moved before the simulation. 


Dependent Storage Elements. The two dependent storage ele- 


\1 12 


m 
| (1) causal conflict 


—sAl > MT Fr—4 | 


\1 m Di2 


i oe 


—A|-— MT FR | + 


(2) algebraic loop, gain 


Di1 m 12 


(3) algebraic loop, gain 


—a |—aAMTF —aA lr 


Fig. 15 Two possible remedies for conflicting /-elements 


ments give rise to a causal conflict, ie., a computational dif- 
ficulty. In the case of the transformer having a constant trans- 
former ratio m one of the J-elements would be transferred to the 
other side and the two of them could be replaced by one eiement. 
As m in the above figure is state-dependent this solution is im- 
practical. Derivative causality is needed, for which 2 solutions 
exist. If it is supposed that differentiation cancels integration, 
both cases give rise to algebraic loops, having inverse gains. 
The loops are easily recognized by following the ‘‘causal arrows’’ 
as explained in the next section. In the numerical implementa- 
tion, however, the differentiation procedure needs an extra unit 
step dealy (see Table 1). Such algebraic loops converge if the 
loop gain < 1. This means that the smallest J-element, cor- 
rected for m?, should have derivative causality. 


Discontinuities. The discontinuous characteristic of a check 


| 
— D+ - wrt 


Fig. 16 Conflicting /- and R-element 


valve or diode forces its causality to be conflicting with the in- 
tegral causality of the inertia element. In such cases a good 
solution can be obtained by replacing both elements by one 
element (J,) having the combined characteristic: 


Ir 
elt 
——4|+—-- 


Fig. 17 Replacement by an input-limited /-element 
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The transmission resistance of the check valve has been 
omitted from the figure, but can be added in the usual way 
without difficulty. Integration within the J, element should be 
executed in such a way, that negative values of e never cause 
negative values of p. It is not acceptable to fulfill this desire by 
clipping negative values at the oudput of a normal integrator. A 
special limited integrator (LMT) as used in THTSIM is needed, 
which zeroes the inpul of the integrator when specified output 
limits are reached, in this case when p <0. 

It is possible to represent J, by the two indicated function 
blocks, because bond graph input to the simulation program mixes 
freely with block diagram input. 

Causal Uncertainties or Algebraic Loops. We consider the fol- 
lowing structure, Fig. 18: 


a ve See 


Fig. 18 Bond graph and equivalent mechanicai network 


Computational causality is dependent on the constraints at both 
ports. The four combinations of constraints are depicted in the 
two causal realizations (4(a) and 4(0)) 


forthcoming. In case 4, 


seem possible, a situation which might be called “causal un- 


certainty.” 


Cc G | R 
2 Ww \ ‘4 


(2) 


T 


—40 41— 


Cc R | R 
a” fee 
| | Gi=1/R1 


pce cereennioomnt G2=1/R2 
Cc RT| G2 
| if 


-+——0 ———lr— -- 


4b 
+—_O—_——— 


Fig. 19 Four combinations of constraints 


A series of ‘causal arrows,’’ —4 , pointing in one direction 
like a one-way road, gives rise to a closed loop of signals in a 
block diagram or signal flow graph. Brown [9] has given an ex- 
tensive analysis of signal loops in unmeshed and meshed bond 
graphs. In Fig. 19, case 1, 2 or 3, no causal one-way roads be- 
tween two algebraic elements exist, though plenty of dynamic 
loops, that include a storage element, are found. 
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By following or plotting (Fig. 20) the causal arrows in case 4, 
however, both realizations (a) and (b) show loops between the 
algebraic elements R; and R2. This is an indication of an implicit 


function or algebraic loop. Such simple, but nevertheless annoy- 


Fig. 20 Algebraic loops in case 4 


ing, algebraic loops can be solved by introducing an iteration 
procedure involving several steps for each time step of the pro- 
gram. In THTSIM a somewhat coarse procedure is followed by 
inserting a block having a delay of one integration step in series 
with one of the two elements, making the iteration time depen- 
dent. To converge, the iteration procedure requires the algebraic 
loop gain to be negative and of absolute value <1. So, in a specific 
case no choice exists between case (a) and (5). 


Conclusion 


It has been shown how block-oriented simulation languages 
The user 
of the program has to augment the bond graph with causal 
strokes. After that causal conflicts and simple algebraic loops 
have to be remedied, as indicated. The THTSIM program with 
bond graph extension is used since 1974 in graduate courses for 
mechanical, electrical and agricultural engineering students and 
for post-graduate courses in modeling and simulation. Students 
have to work with CSMP as well but a strong preference exists 


can easily be extended to accept bond graph blocks. 


for the dedicated minicomputer with its hands-on facilities. In 
research the program is successfully used with large models of 
greenhouse climate control and of pneumatic grain transport. 
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Method of Optimizing Its Convergence 
Properties: 


This paper concerns a generalization of the design of the adaptive observer, and an 
investigation of its convergence properties. By application of hyperstability theory the 
design freedom is enlarged to permit use of time variable adaptive gains. Based on the 
assumption of periodic excitation, and with the aid of Floquet Theory a new convergence 


criterion is defined. By this means the effects of adaptive gains and input frequencies 


upon convergence time are evaluated. 


Structural properties of the adaptive observer 


are compared with other adaptive schemes, in particular with regard to bias due to noise. 


Introduction 


The concept of adaptive identification and observation has 
recently attracted the interest of a number of investigators [1-3].* 
Whether by application of Liapunov or Popov’s hyperstability 
theory, these identifiers are all designed to possess the common 
feature of global stability. Important differences in the various 
designs which have been proposed relate to structure of the 
identifier. In certain cases, as with the adaptive observer, it is 
possible to arrive at a linear representation of the process. As 
will be seen this allows the application of Floquet theory and an 
approach to the convergence problem. 

This paper is concerned primarily with the adaptive observer 
although results are shown to be of more general interest as well. 
In previous work [2-3] the adaptive observer was shown with 
the aid of Liapunov theory to be asymptotically stable in a 
noise free environment, subject to the restriction that the adap- 
tive gains are assigned constant values. By application of hyper- 
stability theory, this result can be generatized so as to remove 
this restriction. As has been shown by Landau [4], this gen- 
eralization in the design is applicable to other identifier con- 
figurations as well. 

Given a stable design for the adaptive process the important 
factors such as rate of convergence and noise response charac- 
teristics must be considered with regard to rate of convergence, 
methods have been proposed for speeding up the adaptive 
process [il]. However no practical method of determining an 


1This work was supported by NSF Grant #37841 and the University of 
Connecticut Computer Center Grant GJ-9. 
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optimal set of adaptive gain parameters and input-signal charac 
An approach to this problem has 
been formulated by recognizing that the adaptive process can 


teristics has been reported. 


in some cases be characterized by a linear periodic differential 
equation—this being a characteristic which can be identified with 
the adaptive observer. In this way, based on Floquet theory [5], 
the convergence rate can be related to the eigenvalues of a dif- 
ference equation, and a meaningful cost function can be de- 
fined [6]. 

Computer studies are reported which demonstrate the de- 
pendence of convergence rate upon values of the adaptive gains 
and input frequencies, and in addition indicate the benefit which 
can be derived from the use of time variable gains. 

The important problem of bias appearing in some of the pa- 
rameter estimates as a result of measurement noise is discussed 
with respect to particular identification schemes. 


1 Deterministic Schemes for Adaptive Identifi- 
cation and Observation 


In this section the basic structures for adaptive identification 
and observation in continuous time which have been. reported 
are discussed with the twofold purpose of comparing, the respec- 
tive design (see Section 5), and placing the adaptive observer in 
perspective. Results concerning the design of discrete-time 
identification and observation will not be discussed in this paper. 
In ail cases to be treated it is assumed that the plant is repre- 
sented by a linear constant ordinary vector differential equation 
of known order, and that the plant is stable, controllable and 
observable. 

In the simplest. identification scheme shown in Fig. 1 it is as- 
sumed that all state variables are measurable. 
essentially a model following configuration in which the param- 


This scheme is 
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Fig.1 Simple identification scheme 


eters of the model are adjusted so as to identify the A and B 
matrices of the plant. The adaptive law relating to this situation 
can be formulated in two different ways, each having certain 
merits and defects which will be discussed in Section 5, 

The scheme in Fig. 2 is similar to Fig. 1 with the important 
distinction that state variables of the plant are assumed to be 
unmeasurable. For this case a stable adaptive identification law 
has been derived [1] by employing state-variable filters [7]. Due 
to the time lag inserted by the state variable filters, present state 
information is unattainable with this configuration. 

The scheme in Fig. 3 represents the adaptive observer [2-3], 
where again the state variable filters are used to enable the de- 
velopment of a stable adaptive law. In this case, however, the 
observer is capable of identifying and observing the plant simul- 
taneously as a consequence of the fact that the filters are not 
placed in series with the plant or observer. The point to be 
stressed is that the location of the state variable filters is a struc- 
tural feature of importance. 

In each of these cases it is found that the design can be pursued 
either by employing the classic Liapunov theory, or Popov’s 
hyperstability theory. The latter which has been used exten- 
sively by Landau [1], e.g., for the configuration in Fig. 2, leads 
to a more general solution to the adaptive problem than has been 
obtained using classic Liapunov methods. In this paper similar 
results are obtained for the adaptive observer depicted in Fig. 3. 


3 Generalized Design of the Adaptive Observer 


The representation of the adaptive observer to be used here 
is taken from [3] wherein the observation-error equation is written 
in a form which proves to be amenable to the application of 
hyperstability theory. 

In the following discussion it is assumed that the plant has a 
single input and a single output. 

Let the nth order plant be defined by 


x = Ax + bu 
y = efx, cT = (100... 9 (1) 


where by observability it can be assumed that A has the form 
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Fig.2 Identification schemes for systems with unmeasurable states 
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The observer is in turn described by 
z= Kr+(k—aén+bu+o 


where 
| J 
K = —k 
0 


is a stability matrix, a, b are estimates of a, b, and } is composed 
of signals added to assure stability of the adaptive process. The 
error equation in turn becomes 


¢é= Ke+ay+6u+o (3) 


wherein e 4 z — x defines the observation errors and « Aa — 4, 
§ Ab — b define the parameter estimation errors. @ consists of 
some auxiliary signals which approach zero asymptotically with 
time [3]. 

By employing the scheme shown in Fig. 3, with state variable 
filters? G(s) defined in the s domain by 


V(s) = G(s)¥(s) 
Q(s) = G(s)U(s) 
it has been shown [3] that an equivalence to (3.3) can be written as 
¢ = Ke + d(v7a + q76) 
€ = c7e (4) 


with e, = €, where e; and €; are the 7th elements of e and e, 
respectively. If r 4 (a?v + 6%q) is considered as the input in 
(4), then the transfer function of (4) becomes 


H(s) 4 e?(sI — K)"d. (5) 


In previous works [2-3], by using Liapunov theory and the 
Kalman-Yakubovich lemma, an adaptive law has been derived 
which guarantees asymptotic stability. However, the adaptive 
gains were specified as constants. This requirement can be re- 
laxed by using hyperstability theory. 

Hyperstability requires a condition on the inputs to the sys- 
tem. Considering the system (4), hyperstability requires that 
the state vector € remains bounded if the input is restricted to a 
subset of all possible inputs. This subset is defined as those inputs 
which satisfy the inequality 


lle(t)ll, for alle >0 (6) 


f r(tje(t)dt < dille(O)ll sup 
6 0 ; ee T 


3An important property of G(s) from a practical standpoint is that lim 
G(s) /s =0. (3s) 7 © 
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Fig. 4 Application of hyperstability theory to the adaptive observer 
in proving the convergence properties [6] 


where 6; = 6;(£(0)) is a positive constant. The system (4) is de- 
fined as hyperstable if, for any r belonging to this set, the in- 
equality 


lle()ll < A(ie(O)il + 61) (7) 


holds for some constant k, and all ¢ > 0. 
If the system (4) satisfies (3.6) and (3.7), and 


lim e(t) = 0 


to 


then the system is termed asymptotically hyperstable. 
With the definitions, we have the following 


Theorem | [8]. A necessary and sufficient condition for a trans- 
fer function H(s) to define an asymptotically hyperstable system is 
that H(s) be strictly positive real (SPR), in the sense that H(s — o) 
is positive real, where o is a sufficiently small positive number. 

It should be noted that this theorem provides the condition 
for satisfying (7). Hence if H(s) is SPR, then (6) can be replaced 
by 


f r(tje(t)dt < 6, for all r > 0. (8) 
0 


The condition of the theorem can be satisfied since it is shown 
in [3] that there exists a nonunique vector d such that H(s) in 
(5) is strictly positive real if K is a stability matrix. It remains 
to be shown that the input r(¢) can be made to fall within the re- 
quired subset. With reference to Fig. 4 which depicts the sys- 
tem (4) together with an adaptive law, this will be accomplished 
by an appropriate choice of W, as indicated. 

It will first be shown that (8) is satisfied; that is, that there 
exists a constant 6 such that 


> 
J r(A)e(A)dA < 6? for all ¢ > 0. 
0 


Using the notation defined in Fig. 4, 


t e t 
J red\ = J ned + f r2€,0X. 
0 0 0 


Considering the first term in (9) 
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t t 
f ned = f aTvedr 
0 0 


If the M matrix is assumed to be nonsingular with M 
(m, M2, ..., Mn), then it follows that 


v(A)e(A)dA = — M-lda(A) 


and by a change of variable (10) becomes 


t a(t) 
f néed\ = — f 
0 (0) 


a) 5 
= f ayme day. 
a® i=] 


0 < & < mz for all i, A. 


a™M—da 


Now let m;~1 be such that 


Then one has 


a(n a 7 a(t) — 
f = am: da; > nf y aida; 
—_ 
(0) a(0) 


=i Lae | 


> _ Y (are) — a0). 


i | 


Since a:(0) can be bounded arbitrarily, it is always possible to 
find a finite 6,2 such that 


hs, 
3 by (att) — a(0)) > — dr. (14) 


t=] 


In the same manner, with N A diag [m, m, 
found so that 


++, Mn], & 5: can be 


k n 
5 by (BA) — Bx(0)) > — 5. 


i=l 


Let 6? A 6,2 + 6,2. Then one concludes that 


t 
J red < & forall t > 0. (15) 


Thus it has been shown that, with the proper choice of d, K, M, 
N, the system defined by 
¢ = Ke + d(v7a + q’6) 
@ — Mvé 
6 — Naé 


is asymptotically hyperstable, and therefore that lim ¢— 0. This 


t7o 


(16) 


is essentially the result obtained in previous works, with the 
added generalization that the adaptive gains { mi, ns} can be 
time dependent quantities, the only stipulation being that they 
be positive. 

As shown in [2], if the input r contains at least (nm + m)/2 
sinusoidal components, where n and m are the nuinber of param- 
eters in A and b, respectively, to be identified, it follows that 
lim @ = lim § = 0. 
to to 

For future reference it will be helpful to rewrite (16) in the 
form 


E = A(E (17) 


with 
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ET A (eT, at, 67) 


ray 
Ae —Mv {0 


—Nq ;0 


As seen in Fig. 3, the signals v, q are obtained by passing the 
input and output of the plant through state variable filters. 
Hence they are explicit functions of time. Therefore (17) char- 
acterizes an unforced linear time-variable system. 

The problem is now to develop a method of optimizing the 
convergence of &(¢) in response to an arbitrary initial unknown 
state &(0). It is noted that v(t), and hence A(t), are dependent 
on the initial state of the plant unless the plant is in the steady 
state. Therefore in order to have a well posed optimization 
problem it is necessary to assume that v(¢), and hence the plant, 
is in the steady state. Failure to make this assumption would re- 
sult in having to find an optimal solution for each initial condition 
of the plant. 

Since very little can be said about the convergence properties 
of a linear time variable system in general, it is expedient to 
make the additional assumption that v(t), q(t), \/(t), N(¢) are 
periodic (AJ, N may be constant matrices), so that A(/) = A(t 
+ T). 
thereby enabling the convergence problem to be cast in a work- 
able form. 


By this means it is possible to apply Floquet theory, 


3 Adaptive Observer—Periodic Case 


As noted above, by assuming that the plant is in the steady- 
state, it is possible to construct A (¢) in (17) so that it is periodical- 
ly time varying. This is advantageous because it then becomes 
possible to characterize the system by a constant matrix. Based 
on the theory of Floquet we have the following 


Theorem II [6]. Given a linear periodic system, 


x(t) = A(t)x(t), ¢ > 0 (18) 


where A(t) = A(t + T), and A(t) satisfies the conditions of con- 
tinuity of solution, then A(t) is globally asymptotically stable, in 
fact exponentially stable, if the eigenvalues of a constant matriz, 


C, lie inside the unit circle, where C is defined as 
CAC,\ra0' (19) 
and 
C,ADRE+1)T+7,kT7+ 7,057 < T 


where B(, ©) is the transition matrix of (18), and k is a non- 
negative integer. 
time t 


Furthermore, with a general choice of sampling 
© 7f+kT.&k = 01,2, ..., the response of (18) is defined 
at these sample points by the solution to the constant difference 
equation 


x(k + 1)7 + 7] = Cx(AT + 7),0<7 <7 (20) 
with the initial condition x(r) being the solution to (18) att = 7. 

By using well known properties of the transition matrix, and 
the assumption of periodicity of A(é), it is possible to show an 
important property of C,;, namely that. [12] 


Cr = P(r, 0)CP(r, 0)-1. (21) 


Thus C; is related to C by a similarity transformation, and the 
eigenvalues of C,; are seen to be invariant with 7. 

By design it is known that {A(c)}, the eigenvalues of C, lie 
inside the unit circle. It is to be shown, with the aid of Theorem 
II and (21), that a meaningful convergence criterion can be 
formulated in terms of a function of {rX(c)}. 

Let the proposed convergence criterion be defined as 


p S max |A;(c)! 


‘ 
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Fig.5 Justification of p (convergence criterion) 


, n, are the eigenvalues of C. To justify 
the selection of (22), consider first the solution to (20) for 7 = 0. 


1, 2, 


where A;, 7 = 


If C and x(0) are known, then x(/) is completely specified at 
t= kT, k = 0,1, ..., and has the form‘ 


x(kT') = = cA. 


as | 


(23) 


Here the vectors of ¢; are determined by the initial condition 
x(0). Since it is desired that a convergence criterion be chosen 
which is satisfactory for any initial condition, it follows that con- 
siderations should be given to the choice of x(0) which causes 


the ¢; associated with max |\;| to dominate the terms in (23). 


Hence (22) is an appropriate measure with 7 = 0. However it 
should be recognized that the behavior of x(t) in the interval 
ee 


there is no a priori knowledge of the behavior of the system over 


T depends upon the initial condition, x(0), and that 


this period except that x(¢) is continuous. It will be shown that 
unless there is a large change in x(t) in this first interval the 
criterion p can be used effectively as a prediction of the con- 
vergence time. Consider the following situation depicted in Fig. 
5 where x(t) is a scalar function. In the first period, the peak 
occurs at ¢ = t; where |x(t,)| > |2(0)| 4.1. By (20) for rT = 0 
a[(k + 1)T] = Cr(kT) (24) 

and for tT = 

a[(R + 1)7 + th) Cax(AT + th) 

= C2(kT + t). (25) 


After N periods 
from t = 0, the ratio of convergence will be C’. That is 


Here C,; = C because x is a scalar. Consider (24). 


a(NT)/x(0) = CY. 26) 


Now to compute the time, &., which will be required for |(t; 
+ &7')| to be less than or equal to |C’x(0)| it is necessary to find 
kT such that 


lajkT + &)| = |C*x(t)| < |C%x(O)|. 


From (27), it follows that 
|CN-k| > |x(t,)/x(0)|. 


Now define the ratio |.xr(t;)/x(0)| A p. The larger the value of p, 
the more severe the degree of variation of x(¢) in the first period. 
The purpose is now to establish a relationship between p and the 
number of periods (&) required for the ratio, |:r(¢.)/2(0)|, to be 
less than or equal to |C’| [Note that C in this case is the eigen- 


value.] We see from (28) that 


k In (1/p) 
_ 


(29) 


N- * Nin(c]° 


‘For simplicity the case of repeated roots is not considered here. The more 
general case can be treated, however, at the expense of greater tedium. 


Transactions of the ASME 





Equation (29) expresses a measure of the accuracy of the con- 
vergence criterion. For high accuracy k/N should be close to 
unity. Consider the case of interest, for which p > 1. Then the 
second term in (29) will be positive. For high accuracy of the 
convergence criterion this term should also be small compared 
to one. The observation is made that p can be relatively large 
due to the logarithmic relationship in (29) without causing k/N 
to be significantly greater than one. 

To generalize this result it must be assumed that if x(0) is 
chosen so that, for 7 = 0, the solution (23) is on. an eigenvector 
corresponding to A;, namely 

x(KT') = ¢;A;*, (30) 
where A; is the eigenvalue having the longest magnitude. This 
then represents the “worst set’’ of initial conditions, in other 
words, the set which causes \; to dominate the response. Since 
it is desired that the optimum design be satisfactory for any 
initial condition, the min p criterion is meaningful for the general 
case. Although (29) has not been generalized, simulation results 
confirm that the criterion can be used successfully in predicting 
convergence. 


4 Optimization of the Adaptive Observer 


Within the constraints of hyperstable asymptotic stability 
there remains some design freedom in the specification of A(f). 
K and d must meet the condition that H(s) 
in (5) be SPR; W(t) and N(é) are specified as diagonal matrices 
with elements which are positive but may be time dependent; 
v(/) and q(/), as output of state-variable filters, must have suf- 
ficient frequency content to guarantee that the eigenvalues of C 
lie inside the unit. circle. 


By way of summary: 


(Failure to meet this requirement can 
result in one or more eigenvalues lying on the unit circle.) A 
complete study of A(é) with respect to all of these variables 
would be necessary for a thorough treatment of the convergence 
properties of the adaptive observer. However, in the work re- 
ported here the study is limited to optimization of the adaptive 
gains, and the amplitude and frequency of the signal components 
entering the plant. 

Intuitively it is obvious that the frequency content of the plant 
input will have an important effect. Also, experience has shown 
that convergence rate is highly sensitive to the selection of the 
adaptive gains. The reason for considering time-varying adap- 
tive gains here is that thereby the design freedom is enlarged. 
This point can be clarified by the following example. 

Consider a simple system as could arise with a second-order 
plant with one unknown parameter. 
(17) is defined by 


Here the error equation 


1 »v 
kn 0 0 
m 0 0 


ku = — 0.25 
ka — 0.30 
v = (0.2) sint 


m = m + mysint 


Letting m: and mz be the design variables, we minimize the con- 
vergence criterion, p, with respect to m, and me, for the cor- 
responding constant system matrix C as described in Section 3. 
To verify the claim of an improved design, the initial points for 
the optimization were taken along the constant gain axis, that is, 
along m. The idea is that if the minimum cost occurs on the 
m, axis, then the time varying part, mz: sin ¢, does not contribute 
to better convergence. However, if m2 turns out to be nonzero, 
then the time-varying part has made a contribution. The com- 
putational results are shown in Fig. 6. It is shown that the mini- 
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Fig.6 Results of gain optimization (note that M;~0 always for mini- 
mum p) — 10 points on ™; axis are the points before optimization; 
After optimization (minimum ») these points are proved to upper 
right hand corner 


mum cost does not fall along the m; axis. It will be seen that the 
contribution of the time varying part is more noticeable as the 
number of adpative gains increases, that is, as the number of 
unknown parameters increases. 

The remainder of this section wiil be concerned with a case ex- 
ample to show the results of optimization with respect to the 
various design parameters. 


Example. The plant is given in the following form: 


-—5 1 Ba 1 ju 


4+ 
—-10 O 
y = [1 0) 
Ze 


(31) is completely controllable and completely observable and is 
represented in the canonical form suitable for an adaptive ob- 
server design. An appropriate adaptive observer is defined as 


z= Kz+k—-—aWJun+bo@ut+wte 
a = [1 Ojz 


a(t) = 


Following the design approach outlined in [3], 
0 r 0 
ar r=-e eee 
“MDw | " [ a™NDiq 


with 


dh: 


D, = and d= 


dy 
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Tatel v ard q are outputs of state variable filters and M, N are constant 
diagonal matrices to be optimized. It should be noted that the 
pair [K, d2] is chosen such that K is a stability matrix and 
e7(s] — K]"ld. is SPR. The adaptive process should be such 
mi me that as ¢ increases, the following results are obtained: 

€ e 
rag on ae ne 15Ad tr ‘ § lim a(t) = 5 lim 6:(¢) = 1 
(pairs are slaved) tro ge 


Varia les to ; 
be opt. mized Optimal values Pmin™ 


: lim a(t) = 10 lim bo(t) = 2 
m2(t)® ’ to tc 


0.025 -—0.19 —0.4 


m(t) = 
m(t) = ne(t) 


and lim z = x. (33) 
(2) These values are normalized based on 7’ = 2II ete 
b) _ - Te \]-1- - 7 
| o ~ es ek ee a ey For asymptotic stability the input u should contain at least two 
distinct frequencies. 
Table 2 Table 1 summarizes the results obtained when p was minimized 
, as a function of the adaptive gains, for the cases in which M, N 
Optimal values of are constant and time varying. Here the input was chosen as 
Assigned Values the variables u = 5.0 (sin ¢ + sin 2.5 ¢). 
1 ) K Ko Pmin™ 


Table 2 summarizes the results obtained when p was minimized 
a. 0.65 as a function of K,, Ke for various values of @;, we with wu = K, 
32. 0 62 sin wt + Kesin wet. Here M and N were assigned a set of nominal 
ae ry values from the top of Table 1. Since the period T was chosen 
70; differently for different frequencies, the value of p has been 
normalized to a standard period. 
= : : The time plots in Fig. 7 are samples of the response of one of the 
44. ‘5 parameter estimation errors corresponding to the data in Table 
44.8 A 1, while those in Fig. 7 correspond to data in Table 2. Com- 
. 35.6 0. parison of p (A max |\;i(c)|) for the various cases shows that the 
16 3. 31. ; 
32 36. 46. ; frequency content is more significant than the adaptive gains in 
16 “hf 45. minimizing the convergence time. However in order that the 
32 7s 29 data not be misleading it should be noted that the system de- 
= fined by (17) is of sixth order. Hence the fact that p is very 
32 A nearly zero signifies that the convergence time could be as much 
(@) These values are normalized based on as 67’, depending upon the initial condition &(0). 
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Fig.7 The piots of parameter identification, b: (see Table 1) 
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5 Effects of Noise 


It is well known that the presence of noise in adaptive proc- 
esses can give rise to bias errors. In the following discussion, 
design factors which can be used to ameliorate this problem will 
be discussed. 

The bias problem has been discussed in [10]. There it is shown 
that, given an adaptive law dit) = ke(t) W(t), if E(d| = 0 in the 
steady state (as is required for stability), then if e and W contain 
correlated noise, the steady state value of $(¢) will contain a bias. 

To see how this result pertains to the present problem consider 
first the identifier in Fig. 1 in which the full state vector of the 
plant is assumed to be measurable. Cases 1 and 2 are two alter- 
native design formulations, the essential differences being in the 
descriptions of the model, and the consequent forms of the error 
equation. 


Case 1 [9] 
Plant x = A,x + b,u 
Model x = Ax + bu + C(% — x) 
Error Equation e = Ce + (A — A,)x + (6 — b,)u 
Adaptive Law a; = — e7Pxi, 
b= — e’Pu 


i = 


Amy scag @ 


Case 2 
Plant x = Apx + b,u 
Model x = Ax + bu 
Error Equation ¢€ = Aye + (A — A,)k + (b6 — b,)u 
Adaptive Law @ = — e’PZi, 1 = 1,2,...,0 
b = — e’Pu 


Tere a; is the ith column of the adjustable A matria of the model, 
and e Ax — x. The distinction between the adapiive iaws in 
these two cases stems from the formulation of the equations for 
the model. 

In Case 1, by forming the model equation as shown it is pos- 
sible to obtain an error equation in which C is an arbitrary 
stability matrix, and the plant state (x) is a multiplier of (A 
— A,). As a consequence it can be shown that the plant state 
variable (z;) appears in the expression for a;. On the other hand, 
in Case 2 the model equation is written so that the error equation 


contains A, in place of C, and the model state (x) is a multiplier 
of (A — A,). This leads to an expression for a; in which the model 
state variable (x) appears. 

The weakness in Case 2 lies in the fact that stability of the 
adaptive process is assured only for a given A, (or at best a con- 
servative estimate of A,). To appreciate the advantage of Case 
2 over Case 1 it is necessary to assume that noise is present in the 
measurement of the plant state x. The terms in a; of Case 1 
then contain correlated noise leading to bias errors in a;. In 
Case 2, however, the noise level transmitted to 7; can be greatly 
reduced, with the result that bias errors in a; can be made smaller 
than in Case 1. An additional observation to be made is that the 
system in Case 2 cannot be represented as a linear system as 
with Case 1. Therefore the optimization technique which has 
been presented cannot be used. 

These conclusions are also in essence applicable to the systems 
in Fig. 2 and Fig. 3 in which measurement of the state vector of 
the plant is not available. Summarily, the identification scheme 
in Fig. 2 is analogous to Case 2 whereas the adaptive observer is 
analogous to Case 1. Thus although the identifier in Fig. 2 is 
superior with respect to bias errors, stability is guaranteed for 
only a small uncertainty in the A matrix, and the system cannot 
function as an observer. The identifier in Fig. 2 also fails to 
qualify as a linear time-varying system, to which the optimization 
technique of this paper can be applied. The adaptive observer 
(Fig. 3) is inferior with respect to bias errors, but functions 
simultaneously as an identifier-observer, and is stable providing 
A is a stability matrix. 


6 Sensitivity Analysis—The Auaptive Observer 


For practical purposes, it is important to know how sensitive 
the optimal adaptive observer is with respect to the parameter 
variations of the unknown plant. The fact that this important 
question does not appear to have been addressed before with 
respect, to the various adaptive schemes which have been dis- 
cussed, is presumably due to the absence of an optimization cri- 
terion. 

Suppose that the parameters of an unknown plant are dif- 
ferent from the nominal values for which an optimal adaptive 











Fig.8 Time plots of parameter identification, a: (see Table 2) 
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Fig. 9 Optimization sensitivity of parameter identification 


observer has been designed. Then the adaptive observer will be 
suboptimal and the speed of convergence may be unacceptable. 
By using the criterion of optimality (p) which is the magnitude 
of the maximum eigenvalue, it is of interest to investigate how 
sensitive this criterion will be with respect to variation of the 
plant parameters. In Fig. 8 are the plots of p with respect to 
individual plant parameters when deviated from their nominal 
values, i.e., the values for which an optimal adaptive observer 
has been designed. It is noted that the variation of an unknown 
plant parameter does not induce any instability of the adaptive 
process unless the plant becomes unstable, in which case the 
system is no longer useful from a practical viewpoint. 

An objective evaluation of these results requires a comparison 
with the sensitivity problem encountered in a fixed (nonadaptive) 
configuration. The purpose here has been to point out an im- 
portant design evaluation criterion which up to now has been 
ignored. 

In order to cope with this problem, it is noted that off-line 
optimization can in theory be carried out without any knowledge 
of the plant parameters, using only a knowledge of A(/) in (17) 
over one period. Specifically, with reference to (17), it is seen 
that a record of v(¢) and q(t) over one period is sufficient to 
permit an off-line optimization to be carried out. 


7 Conclusions 


The adaptive observer, as well as other adaptive schemes for 
identification and observation, must be designed to have good 
convergence properties if they are to be useful engineering con- 
cepts. In this paper a generalization of the design of the adaptive 
observer to permit use of time-variable adaptive gains is derived, 
and a method of optimizing the adaptive process, based on the 
application of Floquet theory, is given. The optimization tech- 
nique is used to study the sensitivity of convergence time to the 
input frequencies and the adaptive gains. Heretofore this has 
not been a practical feasibility due to the absence of an effective 
design criterion. Structural properties of various identifier con- 
figurations are discussed, with emphasis on the effects of noise. 
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The conclusions are that, while the adaptive observer exhibits 
some desirable properties in comparison with other configurations 
in its present form, it is susceptible to bias errors due to noise. 
This problem is viewed as a worthwhile area for further re- 
search. 
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A dynamic interaction model is formulated for the heave-pitch motion of vehicles 
crossing elevated flexible, randomly irregular spans. 
a vehicle passage is modeled using a Bernoulli-Euler beam model and modal analysis 
techniques. Four types of random irregularities characteristic of elevated guideways are 


Span dynamic motion due to 


modeled numerically including vertical span offset, pier misalignment, camber, and 
surface roughness. Analytical expressions for each irregularity power spectral density 
are derived and the relative contributions of irregularities to vehicle excitation are 
examined. The limitations to vehicle passenger comfort levels posed by guideway 
deflection and irregularity are illustrated for personal and rapid transit types of vehicles. 


Introduction 


In the last few years a number of advanced ground transport 
systems which use dedicated guideways have been proposed for 
both moderate speed, urban, and high speed intercity applica- 
tions including automated guideway transit (AGT), advanced 
rail, and tracked levitated air cushion and magnetic vehicle sys- 
tems [1-3].! The potential for implementation of a system de- 
pends to a significant extent upon its capability to provide safe, 
comfortable, and timely service and also upon cost. For applica- 
tion in and around cities a major fraction of the total system cost 
in terms of both economic and environmental impact is related 
directly to the guideway. In many urban areas and also in some 
rural areas, substantial portions of advanced system guideways 
may be elevated to negotiate rights-of-way and to provide 
safety. To minimize visual impact and cost, small cross-section 
guideways with long, graceful spans are desired which may be 
built without stringent construction tolerances with respect to 
alignment and local roughness. The levels to which guideway 
designs and construction techniques may be modified to achieve 
low cost, aesthetic structures are limited by vehicle-guideway 
interactions. These interactions determine the vehicle loads pro- 
duced on the structure and the safety and comfort provided to 
vehicle passengers. 


‘Numbers in brackets designate References at end of paper. 
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In the past five years analyses and performance data for ve- 
hicle-guideway interactions have been developed in a number 
of studies reviewed in [4]. The studies reviewed and more recent 
studies [5-7] have provided methods of determining the vertical 
forces on and deflections of an elevated guideway modeled as a 
smooth Bernoulli-Euler beam due to the passage of a vehicle 
and the accelerations and motion of a vehicle due to guideway 
dynamic motion. Numerical and some analytical methods for 
characterizing the static random irregularities in a guideway such 
as surface roughness, pier misalignment, and camber which result 
from construction and environmental effects have also been de- 
veloped [8-15]. In a few studies [11, 12] the simultaneous ef- 
fects of selected types of irregularity and flexibility on vehicle- 
guideway performance have been evaluated using numerical 
methods to represent irregularities. 

In this paper a vehicle-elevated guideway interaction model is 
developed which includes both the effects of flexibility and ir- 
regularity. Direct numerical and analytical, power spectral 
density techniques of representing random vertical span mis- 
alignment, random pier misalignment, random and periodic 
camber, and random surface roughness are developed. The char- 
acteristics of each type of irregularity are examined in detail. 
The irregularity models are combined with a Bernoulli-Euler 
beam span model and a two-dimensional vehicle model to form 
a coupled vehicle-guideway interaction model. The dynamic 
performance of vehicles characteristic of small and large AGT 
vehicles crossing the flexible, irregular spans are determined to 
illustrate the constraints placed upon span flexibility and ir- 
regularity levels by passenger safety and comfort criteria. 
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Vehicle-Guideway Interaction Model 
The Vehicle-Guideway System 


The general vehicle-guideway system model is shown in Fig. 1. 
It is formulated to determine the vehicle and guideway vertical 
motion which is excited by the full vehicle weight. (Lateral and 
longitudinal motions are excited by only a fraction of vehicle 
weight.) When only vertical interactions are considered, the 
influence of a vehicle or series of vehicles upon a span may be 
represented by a general force distribution f(z, ¢) which is a 
function of time and may vary spatially along the span, but is 
assumed to be constant across the span width. The force dis- 
tribution is determined by vehicle static and dynamic suspension 
forces. 

The guideway profile presented to a vehicle may be expressed 
as: 


y(z, t) = ys(x) + y(z, t) (1) 


where the total profile y,(z, ¢) is assumed to vary only along the 
span and is the sum of a static profile y,(z) which is a function 


only of location and a dynamic profile y(z, ¢) which results from 
traveling vehicle loads. 

For typical advanced guideway systems, span length to width 
ratios are large enough so that individual spans may be con- 
sidered as beams rather than as plates. In addition, the assump- 
tion is made that the guideway consists of a series of simple, 
pinned-end beams resting upon rigid piers. 


Guideway Static Irregularity Representation 


Irregularities occur in a guideway as a result of construction 
practice, settlement, dead weight loads and environmental condi- 
tions. The guideway static irregularity profile, y,(z), may be 
represented as the summation of 2 number of effects including the 
following four effects illustrated in Fig. 2: 


(a) Span vertical offset resulting from differential span 
height. 

(b) Span angular misalignment resulting from differential 
pier height. 

(c) Span camber resulting from dead weight loading, inten- 
tional precamber, or thermal effects. 





Nomenclature 


a guideway span section area 
A surface roughness coefficient 
Aa(t) modal shape function for mth mode 
G. = mean value of camber amplitude 
probability density 
bs vehicle suspension damping coef- 
ficient 
expectation operator 
EI guideway span flexural rigidity 
generalized force distribution acting 
on guideway 
Tn vehicle sprung mass natural fre- 
quency 
= front (rear) vehicle suspension force 
acceleration due to gravity 
span cross section area moment of 
inertia 
vehicle sprung mass moment of 
inertia about sprung mass center 
of mass 


1 
vehicle inertia ratio = I, / 2 Mola? 


vehicle primary to secondary sus- 
pension stiffness ratio 
= vehicle primary suspension stiffness 
vehicle secondary suspension stiffness 
/-T 
vehicle suspension attachment length 
vehicle suspension pad length 
span length 
vehicle mass/span mass 
Mu vehicle unsprung mass 
My vehicle sprung mass 
an integer 
Fourier transform of p(x) 
camber pulse shape 


n 
P(Q) 
p(x) 

R autocorrelation function 


S., Sa, S-, S,, S, = single-sided power spectral densities 
of vertical, angular, camber, sur- 
face roughness and total irreg- 
ularity amplitudes (respectively) 

time 

vehicle speed 

distance along guideway span 

nondimensional distance along guide- 
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way span 2/l, 
guideway dynamic vertical profile 
due to traveling vehicle loads 
= static irregularity profile 
total guideway profile due to vehicle 
loads and irregularity displace- 
ments 
vertical, angular, camber, and rough- 
ness irregularity profiles 
vertical, angular, camber, and rough- 
ness random variables for the ith 
values 
vehicle front (rear) suspension con- 
tact point deflection 
vehicle front (rear) unsprung mass 
displacement 
vehicle front (rear) unsprung mass 
acceleration 
vehicle front (rear) sprung mass dis- 
placement 
vehicle front (rear) sprung mass ac- 
celeration 
nondimensional guideway deflection: 
y/y* 
midspan deflection due to static load 
of vehicle weight at midspan 
= finite length Fourier transform 
sample interval 
impulse function 
: modal shape function for mth mode 
wavelength 
wavenumber 
cutoff wavenumber 
span natural frequency/vehicle 
sprung mass natural frequency 
3.1416 
span material mass density 
variance of vertical, angular, camber, 
and roughness amplitude proba- 
bility densities (respectively) 
= nondimensional time: wit 
= frequency of mth natural mode of 
guideway vibration 
= span damping ratio 
vehicle damping ratio 


y(z, t) 

ys(Z) 

ye(z, t) 

Yr(X), Ya(X), ye(X), yr(Z) = 


Yviy Yaty Yeiy Yri 


Yosy Yor 
Ys, Yir 
Yas, rr 
Y2f, Yor = 


Yas, Yor 


a4, o2, o2, a 
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Fig.1 General vehicle-guideway model 


(d) Surface roughness caused by local surface variation. 


Because these irregularities result from a wide variety of ef- 
fects including construction tolerances and environmental condi- 
tions, they are represented as random irregularities. Measure- 
ments of guideway profiles [13-15] provide guidelines with 
respect to the amplitude probability distributions and amplitude 
spectral densities for various irregularity types. Using this infor- 
mation, models for each individual irregularity may be formu- 
lated assuming that each irregularity amplitude is an independ- 
ent stationary random process. The total irregularity profile, 
y;(x), is generated by summing the four individual irregularity 
profiles: 


YX) = yo() + yo(x) + ye(x) + yr(z) (2) 


The specification of each irregularity is described subsequently. 


Vertical Irregularity Representation. The vertical irregularity 
represents discontinuous alignment of adjacent guideway spans 
and is characterized in terms of the random variable y,; which is 
the height of the ith span of length /, relative to a zero base 
datum. The irregularity is closely related to manufacturing 
uniformity and installation practice. The profile y,(z) due to 
these effects may be generated numerically by selecting from a 
random number generation algorithm [16]? successive values of 
yvi for each span. Then straight line segments of length l, which 
join successive heights, y,i, are used to construct the guideway. 
The properties of the numerically generated irregularity are set 
by selecting the probability distribution, the mean value and 
the variance of the random number set from which the values 
Yvi are selected. 

The power spectral density of the vertical profile, S,, which 
indicates the distribution of irregularity amplitude power as a 
function of wavenumber 2 (Q = 27/\, where \ = wavelength) 
may be derived as shown in Appendix A for a distribution with 
zero mean and variance ¢,? as: 


o,*l, sin? (0.5QI, ) 


SQ) = 
“) T (0.5Q/, )? 


(3) 


The power spectral density is displayed in Fig. 3 and has a low 
wavenumber asymptote which has a value approaching ¢,2/,/7 
as 22 — 0 and then has a number of lobes which are related to 
span length /,, i.e., at Q = 2rn/l, forn = 1, 2,3... S, 


ae) 


is zero 
and is a relative maximum in between. The local maxima are 
contained in an envelope with an asymptotic slope of the form 
1/Q2, i.e., the power contained in large wavenumber irregulari- 
ties decreases with increasing wavenumber. The total mean 


2The algorithm in [16] allows 


the selection of either a Gaussian or uniform 
probability density. 
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Fig. 3 Vertical offset power spectral density 


square amplitude of the irregularity y.(z) is simply ¢,?. 


Span angular misalignment may 
result from differential pier elevation as shown in Fig. 2 and is 


Span Angular Misalignment. 


commonly associated with initial surveying tolerances and long 
term settlement. For a newly constructed guideway in which 
construction tolerances dominate, this type of irregularity has 
been modeled by a random walk type of process in which each 
pier is misaligned with respect to a previous pier by a random 
error. The guideway profile due to this type of misalignment 
may be generated numerically by selecting numbers from a set 
and adding each number to the previous number in random walk 
fashion to generate the ith pier elevation. Then straight line 
segments are used to join each pair of piers and form ya(z). The 


%Simple pier settlement has also been modeled allowing each pier to settle 
relative to a zero base datum [10). 
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properties of the irregularity are established through the dis- 
tribution type, mean value, and variance of the random number 
set. 

The power spectral density of this type of irregularity has 
been derived in Appendix A for a zero mean distribution with 
variance 4? as: 


Oa? sin? (0.522/,) 


Sa(Q) = Tom, (0.521, 2 


(4) 
The power spectral density is plotted in Fig. 4. For small 
wavenumbers, Q/, < 1, the psd, Sa, is represented by an asymp- 
tote, o.2/mQ2l, which increases with decreasing wavenumber (in- 
creasing wavelength). Because the psd theoretically has an in- 
finite value at zero wavenumber, the total power in the spectrum 
is also theoretically infinite.‘ At larger wavenumbers the psd is 
characterized by a series of lobes with zero value at Q/, = 2an 
= n= 1,2,3.... The peaks in these lobes are enclosed in an 
envelope with a slope of the form 1/2; thus the power at larger 
wavenumbers decreases rapidly with increasing wavenumber. 


Camber Irregularity Representation. An irregular guideway 
profile may result from cambering of spans due to temperature 
gradients, variation in prestress, sagging due to creep and long 
term loading, and intentionally designed precamber. For simple 
spans an approximate model of these effects may be represented in 
terms of the function: 


yet) = Yyeilsin rxz/l,| (5) 


where yi; is the camber amplitude of the 7th span and each in- 
dividual span has a sinwave shape. Due to variation in condi- 
tions, the camber amplitude y.; may vary from span to span. 
Random camber is generated by selecting numbers y-i from a 
random set and then synthesizing the irregularity using (5) with 
the successive values of y-i. By selecting the numbers from a 
distribution with a specified mean value 4,, and a specified vari- 
ance, o-2, the effects of variations from a deterministic precamber 
or uniform thermal camber may be represented. 

The camber irregularity psd has been derived in Appendix A 
for the case in which the random variable set has a mean value 
dG. and variance g,? as: 


sin 0.5(Q/, + 7) 
0.5(Ql, + 2) 


, ol, sin 0.5(Ql, — 7) 

S.(Q) = . a 

4a 0.5(Ql, — 7) 
34.2 


oe 


» 5:(Ql, — 2k) (6) 


k=0 


1+ 


The camber pad is plotted in Fig. 5 for the case of d. = 0, i.e., 
For low wavenumbers, Q/, < 1, the 


ysd approaches a constant value of 4¢2/,/m°, while at higher 
PI 


zero periodic component. 


wavenumbers the psd contains a series of characteristic lobes with 
maxima at integer multiples of span length. The envelope of 
these maxima is bounded by an asymptote of the form 1/024; thus 
the power at higher wavenumbers decreases rapidly with in- 
creasing wavenumber. 

If a camber shape is generated which has finite mean value, 
a. ~ 0, then a periodic component of irregularity occurs which is 
represented in (6) by the summation of impulses 6;(Q/; — 2k7) 
which occur at integer values of span length and which have 
powcr magnitude set by a. 


Surface Roughness Representation. Local surface roughness 
may occur along a span resulting from construction techniques 
and the effects of environment and vehicle loading. Based upon 


measurements of at-grade guideways [13-15] the surface rough- 


4A guideway psd is of interest only over a finite range of wavelengths; thus 
the values of the psd at very small wavenumbers Ql, < < 1 are not of interest 
and in practice the spectrum may be truncated. 
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ness ‘is represented directly by a psd of the form: 


A 
Q? + Q2 


S,(Q) = (7) 


where Q; is a cutoff wavenumber and A the spectral roughness 
coefficient. For values of 2 < Q,, the roughness psd is constant 
and equal to A/Q; 


while at wavenumbers 2 > Q,; the psd is of 


the form A/Q? and decreases with increasing wavenumber. For 
simple span elevated guideways the cutoff wavenumber may be 
selected as 27//, to represent local span roughness. 

A numerical roughness profile y,(z) may be generated which has 
the psd of (7) by selecting numbers from a random set with 
specified mean and a variance and passing them through a first 
order digital recursive filter [17] with a sample interval Az to 
provide the psd of (7).5 


Complete Guideway Profile. A complete guideway static pro- 
file y,(x) is generated numerically by adding the profiles of each 
individual irregularity as indicated in (2). Because the processes 
used to generate each irregularity are independent, it can be 


shown that the complete static profile psd S, is [23]: 


} 


5The interval Az must be selected sufficiently small so that the roughness 
profile interval is a fraction of the suspension-guideway contact length. 


Transactions of the ASME 





S, = S, + So + S. + S, 


and the complete static irregularity psd may be expressed 
analytically. The relative importance of each of the individual 
psd’s in high and low wavenumber regions is illustrated in Fig. 
5(a) where nondimensional envelopes of the psd functions are 
plotted. At low wavenumbers all the psd’s reach finite constant 
asymptotic values except for the angular random walk irregulari- 
ty which increases as wavenumber decreases; thus it tends to 
play an important role in small wavenumber (long wavelength) 
irregularities For long wavenumber (short wavelengths) ir- 
regularities surface roughness and the vertical irregularity tend 
to dominate since they decrease as 1/Q222 while the camber and 
angular walk decrease as 1/0224. A dimensional total static ir- 
regularity psd is plotted in Fig. 6(b). This profile psd was gen- 
erated using values for each irregularity which have been 


selected to be typical of those anticipated in a new guideway: 


(1) A vertical irregularity with zero mean and o, = 0.38 cm 

(2) A random walk angular irregularity with zero mean and 
o. = lil em 

(3) A camber irregularity with zero mean and o, = 0.73 cm 

(4) Surface roughness with a roughness coefficient A = 1.5 
xX 107m 


The plot in Fig. 6(6) shows that for these values of irregularity 
the total elevated guideway static psd is similar to the total 
static psd measured on a newly constructed at grade guideway 
[15]. The summation of all the elevated guideway irregularities 
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tends to generate a total psd which may be approximated rea- 
sonably well by a psd of the form A/{? since strong components 
exist in the individual psds at QI, 1 of the form 1/2? and at 
Q/, >> 1 which may be bounded by an envelope with a form 
1/XP. 


Guideway Dynamic Motion Representation 


A number of analytical techniques have been developed to 
determine the vertical dynamic motion of simple beams resting 
upon rigid piers due to traveling vehicle loads [4]. In this paper 
the modal analysis technique is adopted in which the span is 
considered as a Bernoulli-Euler beam and the time and spatial 
varying motion of the span is represented as the sum of the beam 
natural modes of vibration: 


y(z,t) = )) Am(t)>m(2) (9) 


m=l 


where @,,(z) are the modal shape functions determined from the 

boundary conditions and the homogeneous Bernoulli-Euler beam 

equation and A,,(¢) are time varying functions determined from 

solution of the forced beam equation. For simple spans of length, 
l,, the modal shape function is: 

. mre 

¢,(z) = sin —— 


(10) 


and each modal amplitude is determined from the equation: 


a 
di + 2EmWm 


1A m 2 , 
Che + thn = —— [ "f(a, talz)dz (11) 
dt pal, J , 


By solving (11) for each mode with (9) and (10), the dynamic 
motion profile of a span may be determined due to the force dis- 
tribution f(z, ¢) resulting from a vehicle passage. In the most 
general case because the vehicle suspension forces may depend 
upon the motion of several spans and vice versa, sets of span 
equations and vehicle dynamic equations must be solved simul- 
taneously. 


The Vehicle Model 


The vehicle considered in this study is a two-dimensional rigid 
body vehicle with mass m, and inertia J, which is capable of both 
heave and pitch motion as shown in Fig. 1. The vehicle is sup- 
ported by identical front and rear suspensions which consist of 
a secondary stiffness & and damping %, an unsprung mass mu, 
and a primary stiffness &,,. The suspension force is assumed to be 
uniformly distributed on the guideway over a suspension pad 
length, /,, which may be selected to represent rubber tire con- 
tact area, or a levitated vehicle pad length and is equal to the 
force generated in the resultant primary suspension spring ky, 
whose effective displacement is based upon the distance between 
the pad midpoint and the point on the guideway beneath the 
pad midpoint. Four second order linear ordinary differential 
equations are required to describe the vehicle—an equation for 
the vertical motion of the center of mass, an equation for rota- 
tional motion about the center of mass and one equation each 
for motions of the unsprung masses. 
marized in Appendix B. 


These equations are sum- 


Solutions of the Interaction Equations 


The vehicle guideway dynamic interaction equations consist 
of a set of four second order differential equations for each ve- 
hicle and a set of m « & second order differential equations for a 
guideway of k spans with each span represented by m modes. 
Dynamic coupling between the vehicle and guideway equations 
occurs due to the suspension forces acting on the guideway and 
The 


the guideway deflection profile presented to the vehicle. 
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vehicle and guideway equations may be strongly coupled dy- 
namically if vehicle body and suspension dynamics exert a 
significant influence upon guideway defiection or may be weakly 
coupled if dynamic suspension forces are small compared to the 
static force due to vehicle weight, and thus the suspension forces 
are essentially constant and independent of vehicle dynamics. 
The vehicle front and rear suspension forces on the guideway may 
be expressed in terms of the static weight forces and the front 
jjoy and rear jz, sprung mass accelerations and front js and rear 
jar Unsprung mass accelerations by combining the vehicle dy- 
namic equations: 


(diag + Gor) 

2 My My @g g 

21, (ijer — der) 
Mol q? g 


— Meg 1 Mu Mu Yis 


Fyy = 
(12) 


(Yor + der) 


My My @ g 


Mu Mu War 


21» (Hes — Grr) (13) 
Mela? g 

When the sprung and unsprung mass accelerations are a 
significant fraction of one g, then the suspension forces have a 
significant dynamic component and strong dynamic coupling 
occurs between the vehicle and guideway. However, a primary 
goal in advanced system design is to achieve passenger compart- 
ment accelerations in the range of 0.059 or less. For vehicles which 
meet this goal jjos/g ~ tjer/g ~ 0.05 and for which the unsprung 
mass inertia forees are small, i.e., (muijis)/(mog) ~ (muijor)/ (mg) 
< 0.05, the front and rear suspension forces on the guideway 
deviate from constant values of 


Mau 


Fyy = Fy aed MG 1 + (14) 


My 


by less than 10 percent. When (14) provides a good approxima- 
tion to the vehicle suspension forces, the guideway deflections 
may be computed independently of vehicle dynamics and weak 
dynamic coupling exists between the vehicle and the guideway. 

When the vehicle-guideway equations are strongly coupled 
and suspension dynamic forces are a significant fraction of the 
weight forces, then the four equations for each vehicle and the 
m « k guideway equations must be solved simultaneously using 
techniques such as direct numerical integration. In this case the 
static guideway irregularity profile y,(2) is generated numerically 
and added to the dynamic profile y(z, £) computed at each in- 
stant of time. From the total profile y,(z, ¢) the vehicle motion, 
the suspension force acting on the guideway, and, in turn, a new 
guideway dynamic profile are computed. Direct numerical in- 
tergration yields the guideway and vehicles motions as a function 
of time and has been used in [11, 12, 18] to determine vehicle- 
guideway system performance. 

When vehicle suspension forces may be approximated as con- 
stants given by (14), then (11) may be solved independently of 
vehicle dynamics to yield a periodic dynamic guideway deflection 
profile which may be added to the static irregularity profile to 
yield a total guideway profile y,(z, ¢). Vehicle time history re- 
sponse to this total profile may then be determined by integration 
of the vehicle equations. The use of the partially coupled model 
may lead to considerable savings in computation time compared 
to the fully coupled case as discussed in [7, 18, 19]. 

The vehicle-guideway motion time histories generated by 
either strongly or weakly coupled cases provide data to deter- 
mine guideway stresses and deflections and vehicle motion and 
accelerations. In cases in which guideways with random ir- 
regularities are considered, the time histories may be analyzed 
to determine statistical properties of the time histories. An area 
of particular importance, for example, is that of passenger com- 
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fort which, although difficult to define [19], has often been assessed 
in terms of rms acceleration, rms acceleration in given frequency 
bands, or acceleration power spectral density. The computation 
of these statistical quantities and the length of time history re- 
quired for accurate statistical data are discussed in [17]. 

For systems in which the partially coupled vehicle guideway 
model is valid and in which a linear vehicle model is used, the 
rms acceleration, in frequency bands, acceleration power spectral 
density, and total rms acceleration may be computed in a direct 
and computationally very efficient manner: 


(1) From the solution of (11) to constant traveling forces, a 
periodic guideway profile is obtained which with classical transfer 
function analysis methods yields vehicle acceleration’ due to 
guideway deflection at each frequency. 

(2) The vehicle acceleration psd due to irregularities is deter- 
mined directly from the guideway psd given in (8) and classical 
transfer function analysis methods [17]. 

(3) The total vehicle rms acceleration due to both irregulari- 
ties and deflection may be determined in any frequency band by 
combining the two independent contributions. 


This direct method fer computing rms vehicle accelerations 
yields data in a fraction of the time [23] required by methods in 
which the time history is generated first.8 For many cases of 
interest in transportation system design, this method of compu- 
tation may prove to be useful, at least in preliminary design 
stages when general vehicle guideway parameters are bing 
selected to provide a system which meets a specified passenger 
comfort criterion. However, if either a nonlinear vehicle model is 
used or ‘a fully coupled vehicle-guideway model is required then 
integration methods to determine time histories and subsequent 
analysis to determine statistical properties is required. In the 
following section vehicle performance on a flexible, irregular 
guideway is determined using the direct method to compute rms 
acceleration in one third octave bands so comparison of vehicle 
ride quality with the ISO criteria may be made conveniently. 


Vehicle Performance on a Flexible Irregular 
Guideway 


To illustrate the influence of irregularities and guideway de- 
flections upon vehicle performance a number of parametric studies 
have been conducted. The two vehicles summarized in Tabte 1 
with 1 Hz suspension natural frequencies have been considered— 
a medium size PRT? vehicle running nominally at 26.8 m/s anda 
larger transit vehicle running nominally at 40.2 m/s. 

Elevated spans which contain the levels of irregularity cited in 
Table 1 were designed for these vehicles so that at the nominal 
operating conditions, the vehicles meet a one hour extended ISO 
passenger comfort specification {19, 20] in which the rms ac- 
celeration in one third octave bands is limited to specific values. 
The span structural properties are summarized in Table 1. 

The results of the parametric studies for the two vehicles 
running on the 15.2 and 30.5 m length spans are summarized in 
Figs. 7-10. In each of these figures the rms accelerations meas- 
ured at the front or rear of the vehicles!’ are plotted in one third 
a higher 


octave bands for (1) the baseline design case, (2) 


Sif a partially coupled vehicle-guideway model is invalid or a nonlinear 
vehicle model is used, then direct integration to determine time histories and 
subsequent analysis of the time history is required to determine statistical 
information. 


™While acceleration is cited, any vehicle or guideway response variable may 
be dete-miined. 


8The time savings depends upon the particular system simulated. For the 


systems presented in this paper it is on the order of a factor of five. 
*Personal Rapid Transit vehicle, a class of small AGT vehicles. 


The front and rear accelerations are similar in magnitude in each frequency 
band; therefore, only the acceleration which is a maximum is presented. 
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speed case, (3) a case in which the amplitudes of each irregulari- 
ty are doubled, and (4) a case in which the vehicle with a 0.75 
Hz suspension natural frequency is operated over the guideway 
with double irregularity amplitude. In each figure the triangles 
indicate the acceleration due to both deflection and irregularities 
while the circles represent acceleration due to only deflection. 
The coincidence of a circle and triangle indicates that within the 
designated frequency band irregularities generate negligible ve- 
hicle acceleration in comparison to deflections. 

Fig. 7(a) shows that the PRT vehicle design on a 15.2 m span 
is primarily limited by span deflection with the maximum ac- 
celeration contribution occurring at the vehicle crossing frequency 
v/ly = 1.75 Hz. At this frequency irregularities do not sig- 
nificantly contribute to acceleration in comparison to the de- 
flections. In higher and lower freqeuncy bands deflection has a 
much less significant contribution to acceleration and irregulari- 


Table 1 Summary of vehicle-guideway parameters 
Vehicle Parameters 
No. 1 


r: m/s 26.8 40 
la: m 3.05 12.3 
M, + M.: kg 2180 16,000 
fa: Hz 1.0 1.0 
My/ My! 0.25 0.15 
Kar/kp: 10.0 10.9 
Iy: 1.0 1.0 
ee 0.25 0.25 


Parameter 


Span Structural Properties 
for vehicle No. 1 for vehicle No. 2 
5.2 30. : 30.5 
.71 X 108 53.6 & 108 
3 8! Q, 1510 
: 0.025 ; : 0.025 
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Be 


ties dominate the contribution to acceleration. Fig. 7(b) shows 
that if the vehicle is run at a 25 percent increase in speed, the 
maximum acceleration occurs at the new crossing frequency of 
2.2 Hz and is increased in magnitude to exceed the ISO speci- 
fication. The accelerations due to irregularities in Fig. 7(b) also 
increase due to increased speed and exceed the ISO specification 
at 8 Hz. 

Fig. 7(c) shows the influence of doubling the irregularity 
At the vehicle 
5-10 Hz frequency range the 


amplitude on the baseline acceleration response. 
crossing frequency and in the 
acceleration levels are increased and exceed the specification; 
thus even if the guideway rigidity were increased, the irregulari- 
ties would provide a limit to passenger comfort. 

The influence of doubling guideway irregularities and simul- 
taneously reducing vehicle suspension natural frequency to 0.75 
Hz is shown in Fig. 7(d). With the lowering of the suspension 
natural frequency, the vehicle is capable of running on a guide- 
way with double irregularity amplitudes and meeting the ISO 
specification. For this design case an improvement in vehicle 
suspension characteristics allows a two-fold increase in guideway 
irregularity levels while maintaining passenger comfort. 

Fig. 8 summarizes data for the PRT vehicle No. 1 operating 
on 30.5 m spans. This figure shows that with the longer span the 
maximum deflection limit occurs at the vehicle span crossing 
frequency cf 0.88 Hz and also shows that at this frequency ir- 
regularity contributions to the acceleration are an order of mag- 
nitude less than the deflection. At other frequencies the con- 
tribution of irregularities to the acceleration is at least a factor 
of two less than the ISO specification. Increasing the vehicle 
speed and/or doubling the irregularities in this long span length 
case is possible without exceeding the ISO comfort specification. 
Thus, this longer span length case may have irregularities of 
double amplitude in comparison to the short span length case 
and meet passenger comfort. Since the vehicle-span crossing 
frequency is low, 0.88 Hz, the frequencies in which deflection 
and irregularities make primary contributions to acceleration 
are separated and for an ISO type of specification the constraints 
posed by the two factors may be considered essentially inde- 
pendently. 
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Fig. 7 PRT vehicle operating on a guideway with 15.2 m spans 
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Data for the higher speed rapid transit vehicle crossing 15.2 m 
spans are summarized in Fig. 9. For this vehicle the maximum 
acceleration occurs at 2.6 Hz, which is the vehicle-span-cross- 
ing frequency. In this frequency range significant contributions 
are made to acceleration by both irregularities and deflec- 
tions; thus even if the guideway were made rigid, acceleration 
levels near the passenger comfort limit would occur due to ir- 
reg ilarities. Increasing the vehicle speed by 33 percent raises 
the crossing frequency to 3.4 Hz and yields accelerations which 
exceed the ISO criteria. Doubling of guideway irregularities also 
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leads to increased accelerations in the 1.5-10 Hz range which 
also exceed the ISO limit, while reducing the vehicle suspension 
natural frequency to 0.75 Hz allows the vehicle to meet the 
ISO specification with double amplitude irregularities. 
Additional data are summarized in Fig. 10 for the higher speed 
vehicle crossing 30.5 m spans. These data show that the maxi- 
mum acceleration occurs at the vehicle span crossing frequency of 
1.3 Hz. On this longer span the deflection dominates the irregu- 
larity contribution at the crossing frequency. For this ca:e either 
increasing the speed or doubling the irregularity level leads to con- 
ditions which exceed the passenger comfort specification. Re- 
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Fig.8 PRT vehicle operating on a guideway with 30.5 m spans 
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Fig. 10 Rapid transit vehicle operating on a guideway with 30.5 m spans 


ducing the vehicle suspension natural frequeney to 0.75 Hz 
allows the vehicle to meet the! ISO specification while operating 
on a guideway with double the levels of irregularity. 


Summary and Conclusions 


In this paper a model for vehicle-elevated guideway inter- 
actions has been developed. The model uses classical modal 
analysis techniques to represent guideway dynamic vibration. 
Numerical techniques are used to generate a guideway irregulari- 
ty sample function which is synthesized from four types of ran- 
dom irregularity including vertical offset, pier misalignment, 
camber, and surface roughness. Analytical expressions for the 
power spectral density of each irregularity and of the complete 
guideway have been derived. These derivatives have shown the 
frequency bands in which each type of irregularity has important 
contributions and have shown that for typical guideways antici- 
pated, pier misalignment dominates long wavelength irregulari- 
ties, while vertical offset and surface roughness dominate short 
wavelength irregularities. 

The influence of guideway irregularities and deflections on 
vehicle performance has been illustrated for moderate and 
Data for the smaller 26.8 m/s vehicle 
crossing 15.2 m and 30.5 m spans showed that when acceleration 


was analyzed with respect to an ISO type of criterion, the deflec- 


higher speed vehicles. 


tion generated acceleration components were primarily located 


at the vehicle crossing frequency v//, of 1.75 Hz for the 15.2 m 
span and 0.88 Hz for 30.5 m span while irregularity generated 
acclerations were of importance primarily above the 2.5 Hz 
range. For these cases the irregularities and deflection generated 
accelerations occur in separate frequency bands; thus for an ISO 
type of specification the deflection and irregularities may be 
considered essentially independently. For this vehicle it was 
shown that irregularity levels could be increased by a factor 
of two when the vehicle supsension natural frequency was de- 
creased from 1.0 to 0.75 Hz and the ISO criteria would still 
be satisfied. 

For the higher speed vehicle it was also found that maximum 
acceleration levels occurred at the vehicle crossing frequency 
which was 2.6 Hz for the 15.2 m span and 1.3 Hz for the 30.5 
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span. At the 2.6 Hz crossing frequency both irregularities and 
deflection made significant contributions to acceleration; thus a 
direct. tradeoff exists in span design between irregularity and 
deflection and the two factors must be considered simultaneously 
for an ISO type of design criterion. For this higher speed vehicle 
it was also shown that the level of guideway irregularity could 
be doubled when the vehicle suspension natural frequency is 
reduced from 1.0 to 0.75 Hz and an ISO type of criterion would 
be satisfied. 

In summary, this paper has presented models to evaluate the 
interactions between vehicles and flexible irregular guideways 
and has shown cases in which principally deflections, principally 
irregularities, and both deflections and irregularities limit guide- 
way design. 
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Analysis 


APPENDIX A 
Derivation of Irregularity Power Spectral Densities 


Vertical Offset. The vertical offset psd may be determined di- 
rectly from its autocorrelation function R(x) which has been de- 
rived in [20] as: 


R(z) 


R(z) (15) 


where: 


& = z/l,; = nondimensional distance along span of length /, 
The autocorrelation and psd are a Fourier transform pair and 
thus the vertical psd S,(Q) may be expressed as: 


SQ) = 3 i R(zx)e~i@=dz (16) 
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which yields: 


o,7l, sin? (0.5Q/,) 


S 2 = 
(Q) = (0.5Q1,)2 


(17) 


Angular Irregularity. The angular random walk irregularity is 
simply the integral of a vertical offset type irregularity. The 


profile generated by a random walk yz is: 


1 z 
Ya(r) = i, f Y»(x)dx 
“ 0 


Thus, the angular psd is related to a vertical psd by: 


(18) 


S,(Q) = S,(Q) 


and the angular psd S,({2) is: 


og. sin? (0.521 Is) 


So(Q) = Tom, 501.7 


(20) 

Camber irregularity. A camber irregularity which consists of 
a series of half sine functions with spacing /, and random ampli- 
tude is shown in Fig. 2. To determine the psd of camber with a 
general shape, it is convenient to write the camber profile as a 
series of centered pulses: 


@ 


» Genp(x — nl.) 


n= —2O 


Ye(Z) = 


-where 
Ye = camber profile 
Gen amplitude random number 
nm = an integer 
and where for sin camber 


Tr 


1 for — 1,/2 < x < l,/2 


p(x) = 


p(x) = 0 for x elsewhere (22) 


The Fourier transform P(Q) of the function p(x) is defined as: 


P(Q) = f- p(x)e~i2zdx (23) 


The psd may be defined using the expectation operator E{ } as: 


8.(2) = 1 lim | Z 


qe == E 


)|? 
OL (24) 


where Zz is the finite length Fourier transform defined over the 
sample length 2L by: 


L 
Zi(Q) = f y-(x)e~i2zdz 


L 


(25) 
Using (21) and (23) in (25) and following a development similar 
to that of [21]: 


N 


= PQ) P) aeneini’s 


n=—N 


Z1(Q) (26) 


The spectral density may be derived by using (26) in (24) with 
the result: 


1 N 
- lim EZ IP(Q)|2 | a ene M's 


T N+0 


2/2) « 1 , 
S(Q) = QN + Dh (27) 


n=—N 
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where the sample length 2L is given by: 
2L = (2N + 1)l, 
Upon expanding and noting that: 
Efaciac;} = 0255; + 342 (28) 
where 


variance of camber amplitude probability density func- 
tion (p.d.f) 


mean value of camber (p.d.f) 
= 1 wheni = j 
0 when 0 ¥ 7 


ce ¢c 


2 342 
= = |p + — IPO)? 
tl, al 


N 


dy eints2 (29) 


n= —N 
' 


li 1 
im = 
veo (2N + 1) 


which after considerable manipulation [22] can be reduced to: 


3a2 = 
i, (POMP YL (Ql, — 2ke) 


k=0 


o2 
S(Q) = |P(Q)|? + (30) 
tl, 


where 6; is the Dirac delta function. 
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APPENDIX B 
Summary of Vehicle Equations 


The equations of motion for the vehicle shown in Fig. 1 may 
be derived in nondimensional form as: 


u (O%y 4, OXe\ 5 MM (as ,@ 
; dr? dr? Q, dr 


\ 
M 2§.M dYiy dYy, 
. rete all 4 
+ 53 (War + Yn) = 5 (= = 


" 
+ AY. (Yas + Yir) 


(31) 
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An Efficient Simulation Method for Distributed 
rs. sive. # Lumped Fluid Networks 


Assistant Professor 
of Mechanical Engineering. 


A method is presented for the simulation of fluid networks consisting of uniform dis- 
D. N. WORMLEY tributed elements and lumped dynamic, nonlinear elements. The uniform transmission 
ii elements may be lossy and dispersive. General relationships are derived for their 
Associate Professor terminations in junctions with other elements and/or with dynamic, nonlinear lumped 
of Mechanical Engineering. elements. The basic computer simulation method is efficient in terms of computation 
Massachusetts Institute of Technology, time and core storage requirements in comparison to direct finite difference methods 
Cambridge, Mass. and may be implemented on a minicomputer. Simulation results are compared with 
experimental data for a pneumatic transmission line terminated with a nonlinear 
resistance and for a pneumatic transmission network consisting of three lines of incom- 
mensurate lengths. 


Introduction variety of networks. It has been developed to be computationally 
efficient in terms of both time and memory requirements in 
comparison to the direct numerical simulation methods de- 
scribed in [4-8] so that it is feasible to be implemented on a 
minicomputer. 

The principal elements involved in the simulation are de- 
scribed in following sections while evaluation of the method by 
comparison with experimental data for a multiple line fluid net- 
work is described in a subsequent section. 


The representation and analysis of fluid networks has been of 
continuous interest for a number of years in a wide variety of 
areas including large scale hydro distribution systems, bivlogical 
systems, heating and ventilation systems, power plant exhaust 
gas systems, and fluidic signal processing and control systems. 
Considerable effort has been devoted to the study of complete 
systems and of individual system elements. Studies have been 
reviewed in [1]! which have been devoted to the characterization 
of fluid transmission lines. These studies have been extended to . . 
work which has also been devoted to single line characterization Simulation Development 
(2, 3]. Complementing the studies of individual lines have been 
papers by Tsao [4, 5, 6], Streeter [7], Schaedel [8], Foster [9], 
Franke, et al. [10], and Wright and Mufti [11] which have dealt 
with the modeling and simulation and, in some cases, experi- 
mental study of fluid networks consisting of multiple lines 
coupled to linear and nonlinear terminations. 

In this paper a unified method is developed for the simulation 
of fluid networks which may consist of multiple transmission 
lines with nonlinear, dynamic terminations. The simulation 
method is developed to accommodate a variety of one-dimen- 
sional, small amplitude, acoustic uniform transmission line models 
varying from simple lossless line models to lossy, dispersive 
models and general nonlinear dynamic terminations. The simula- 
tion yields the pressures and flows in the network as functions 
of time. The simulation method may be directly applied to a 


The simulation is developed for fluid networks of the general 
type shown in Fig. 1 in which the network can be partitioned into 
distributed and lumped subsystems. The simulation is de- 


near Dynamic 

Termination 
nped 

Parameter 


>ubsystem) 


'Numbers in brackets designate References at end of paper. 
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veloped by first treating the distributed elements and then con- 
sidering the lumped elements as terminations. 


Distributed Elements 


The distributed elements are modeled as linear one-dimen- 
sional uniform area, rigid transmission lines with through flow 
mean velocity, small compared to acoustic velocity in which the 
= Z.4q 


where q is flow and Z, is the line surge impedance.? When @ is 


variables are pressure p and ‘‘flow-equivalent pressure’ 


defined as positive into a distributed element at both ends, then 
the values of p and g at end 1 may be expressed in terms of the 
values at end 2 as [1]: 


—sinh T 


) eosh T 
” = (1) 


n sinh —cosh T G2 


For convenience in computation, wavescattering variables 
may be introduced where in general: 
1 +1 Pi 
v5 1 


Ui 
—1 qi 
By combining (1) and (2) and rearranging 
0 eT Uy é 
(3) 
eT 0 U2 
The notation used in (1)-(3) in which the flow variable @ is 
defined as positive into the transmission line at both ends is not 
conventional [1]; however, with this notation terminations at 
each end of the line may be treated identically and general 
termination equations may be derived in a convenient form. 
The wave variable form of the transmission line equations in 
(3) is well suited for digital computation; however, additional 
relationships are needed to transform the variables to p and q 


3In this development, the surge impedance is considered constant; however, 
the extension to a dynamic surge impedance is straightforward. 


variables at terminations and junctions. The two causal forms 


of interest may be derived for: 
(1) Pressure Causality as: 
2 —1 
1 —1 
Flow Causality as: 
Ui qi a 
(5) 
Pi vi 

The relationships for the transmission ‘ine element in con- 
venient computational form are summarized in Fig. 2 where + 
signs are required for a g input and — signs are required for a p 
input. The block diagram is symmetric so that terminations at 
each end may be treated identically. 

Digital simulation of a line element requires the time domain 
solution of the operation e~l + u; and requires solution of the 
algebraic termination relationships. These areas are discussed 
in the following sections, 


Digital Simulation of e-' 


Digital simulation requires the convolution of the inverse 
Laplace transform of e~! with the functions ui = (2p; vi), Or 
ui = (29: + uv); thus the convolution weighting function A(t) is 


required: 
(6) 


The convolution weighting function h(t) depends only on line 
geometry and fluid properties and is independent of the line 
is useful to 


h(t) = £ "e-T] 


terminations. For convenience in calculation it 
factor e~! into a pure delay component and an additional con- 


volution component: 


lr = Tas + I'(s) (7) 


so that the output v; due to an input u; for the line segment join- 
ing nodes 7 and / is: 





Nomenclature 


frequency 


lowest frequency of 


matrix defining dynamics of a 
linear terminator 

modified A matrix 

area of terminating orifice 


response 


first propagation constant 

input matrix for linear dynamic 
termination 

modified input matrix 

second propagation constant 

column k of input matrix 

output matrix for linear dy- 
namic termination 

flow coefficient for quadratic number 
orifice 

discharge coefficient for quad- 
ratic orifice 

speed of sound tion 

Ci row / of output matrix 

D = direct transfer matrix for linear flow 

dynamic termination 

D' reduced direct transfer matrix 

element at row / column & of D 

Ey truncation error 


complementary error function 
F = normalized frequency 
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incremental change 


impulse response of finite por- 
tion of semi-infinite line 
subscript denoting 7th portion 
distributed network 
dummy index 
dummy index 
length of line 
dummy index 
number of convolution steps 
reduced acoustic Reynold’s 


number of time delay steps 
number of lines in a network 
number of time steps in simula- 


pressure at a junction 


flow equivalent pressure: 

radius of line 

input vector for linear dynamic 
termination 

modified input vector 

Laplace variable 


time to process boundary con- 
interest ditions 
in impulse time to simulate transient 
delay time for a fluid line 
time to perform convolution 
nominal delay time: C.o/L 
variable moving into line 
coming out of 


wave variable 


variable line 

variable state vector for linear dynamic 
termination 

variable surge admittance of a line 

source admittance 

output vector for linear dy- 
namic termination 

surge impedance 

source impedance 

propagation operator 

modified propagation operator 

time step 

mean kinematic viscosity 

Z4 dummy time variable 

subscript denoting one end of a 
line 

subscript denoting other end of 
a line 

inverse Laplace operator 
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v(t) = rj u(t — Ta — r)ha(r)dr 
0 


ha(t) =< “Hew T] (9) 


A numerical solution of (8) may be written in terms of the time 
increment At as: 


© (i+DAt 
vi(kAt) = b> ui(kAt — Ta — jdt) { ha(r)dr (10) 
J jat 


7-0 
When the time increment At is selected as: 


At = Ta/n (11) 


where n is an integer and where k and n are used as subscripts, 
equation (10) may be written as: 


v(kAt) = )) ws((k — n — j)At)H; 


j=0 


(i+Dat 
f ha(r)dr 
gAt 


The value of H; depends upon the line model used, i.e., the 
form of I adopted. In this paper a form suitable for circular, 
constant area acoustic level, pneumatic transmission lines which 
includes the effects of losses and dispersion derived by Brown 
[12] is used. In this case 


Bi ro 
-— Ai 2T2 
hats) = @¥ eto| 4 a | 


cor? 
alt 


H; = (13) 


A; = 1.48 


fof 


B, = 1.08 L 


In simulation the infinite summation of (10) is truncated at 
j = mso that an approximation is introduced into the computa- 
tion. The accuracy of the computation depends upon the selec- 
tion of m which in effect selects the time inverval over which the 
convolution is performed when At is specified. The error from 
the truncation may be expressed for the Brown model approxi- 
mately as [13]: 


f er f ae eee 
i SS ~ ere [2a 5, | a) 
f ha(r)dr wa 
0 


A more conveniext form for the error estimate valid for input 








Fig.2 Block diagram of a line segment 
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signals whose lowest frequency of interest is f; has been derived 
in {13] as: 
m 1 | 0.040AiF 


(16) 


0.718 
= — ell .5F0-486/v0-514— 5, /N] 
n F E;N 


F 4T afi 
E; desired truncation error 


To achieve a desired truncation error of 0.1 percent, the value 
of m/n required is plotted in Fig. 3 as a function of the Reynolds 
number N and frequency of interest F. 


Terminations and Lumped Element Representa- 
tion 


Static Linear Terminations. Two of the general types of static 
linear terminations for line segments are illustrated in Fig. 4, 
where a line is terminated with a series impedance Z, and pres- 
sure input P, and where a line is terminated with a shunt ad- 
mittance Y, and flow input Q,. In these two cases algebraic loops 
occur in the termination block diagrams which must be solved 
explicitly to achieve a stable solution method. For these two 
cases the solutions in desired causal form may be derived as: 


for series impedance case: 


2p. — [1 — Z./Z ui 
1+ Z./Z. 
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Fig. 3 Convolution time as a function of frequency and Reynolds 
number 





b.) Flow Source 


Fig.4 Pressure and flow boundary conditions 
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for the shunt admittance case: 


2Z ee + (1 — ¥eZev 


1+ Y,2Z. 


uN= (18) 
By using these equations directly the problem of the instability 
posed by the algebraic loop is avoided. 


Dynamic Linear Terminations. When a line element is ter- 
minated by a linear dynamic lumped element system, either of 
two causalities must be employed in coupling the lumped element 
equations with the distributed line segment. If r;, represents the 
output from the ith distributed element port and y; is the input 
to the ith distributed element port, then when r; is a pressure 
(pi), ys is a flow equivalent pressure (g;) and vice versa. The 
lumped element subsystem may be represented by a set of equa- 
tions: 


x = Ax + Br (19) 


y=G+Dr (20) 


where: 
x = system state vector 
A, B = system matrices 
C, D = output matrices 
y = output vector with y; the line input term 
rt = input vector with r; the line output term 


The termination of a line with a general dynamic system 
represented by (19) and (20) is illustrated in Fig. 5. 

Dy, the element in the D matrix relating the line outputr; 
to the line input y;, has been separated from the D matrix leaving 
D' (for which Dy.' = 0). If Dy is nonzero, then an algebraic 
loop is formed which must be eliminated. For the case where y; 
is @ pressure and r; 4 flow, the loop may be solved as: 


Ci 4 D, 
—< 9 é 
1 — Du 1 — Du : 


v= . (21) 


r= yt % (22) 
where r' is the vector r with r, replaced by v; and C;, D; are the 
rows of C and D which correspond to y:. 


If y: is a flow equivalent pressure and rz is pressure, the only 
modifications required are to replace + v by —v in (22) and in 


' 
r 


Finally it is noted that the coupling of the line with the dy- 
namic system may be treated in explicit form for simulation and 
it is convenient to define: 


B,C; 
1— Dn 


B,D; 
B* = es 
e+ 1 — Du 


and rewrite the dynamic system equations as: 


At =A+ (23) 


(24) 


x= A*x + B*r' 





Fig.5 A linear dynamic termination 
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where the modified input r' is a convenient computational form 
and contains +v for pressure and —v for flow causality in 
place of rz. 


Static Nonlinear Terminations. Static nonlinear terminations 
in which a general functional relationship exists between pressure 
and flow occur in a number of cases of interest including fan char- 
acteristics and orifice characteristics. For these cases whether 
the relationship is expressed in tabular or analytical form, an ef- 
fective nonlinear algebraic loop exists which must be solved. If 
the nonlinear function is sufficiently simple, for example quad- 
ratic in nature: 


Pe all Pi = 9: |9:| (26) 


such as is common for orifice flow at low Mach numbers, then an 
analytical solution with appropriate causality is possible: 


1 a ae 
Reel aa 2c E = Vi + 4c|lp, — al |e (pe — %] (27) 


where p, is the pressure on the side of the orifice away from the 
fluid line. Principle values have been assigned to be consistent 
with physical constraints [13]. 

For more complex nonlinear functions or functions in tabular 
form the implicit algebraic loop may be solved iteratively at each 
time step in the solution technique. Such procedures have been 
employed in [13]. 


Junctions and Branching 


Conditions in which several lines are joined at a common 
point with no other terminations may be analyzed directly by 
applying continuity and compatibility relations. If at a junction 
the pressure is assumed to be uniform and the volume flows are 
assumed to sum, then the block diagram shown in Fig. 6 may 
be used to represent a junction of n; lines. At a junction, causali- 
ty requires that one line have pressure as an input to the junction 
and all other lines have flow as an input. The block diagram il- 
lustrates that for an n; line junction, a set of n; algebraic loops 
exist which must be solved to yield each v; component in terms of 
u; components and each line surge admittance Y,; = 1/Z,. 
From continuity: 


" 
Yan = — )) a¥u 
i—2 


and noting the junction pressure as P;: 
Ppeutin 
and 
Gg = Py —% 


Thus, the pressure is: 


Fig.6 Bieck diagram of a junction structure 
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pa v5 Ye; 
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P,; = ——— 


“s 
yy. 


it~] 


and each wu; is: 


uy = 2Py — 4% (32) 
Equations (31) and (32) provide the computational relations 
required to characterize a junction of n; branches. 


Overall Simulation 


The general simulation method is shown in Fig. 7 and consists of 
three types of calculations: 

(1) The static functional calculations represented by termina- 
tions and junction equations 

(2) The numerical multiplication and summation represented 
by the distributed elements 

(3) The numerical integration of the linear differential equa- 
tions 


These three basic steps have been implemented in a digital 
computer program. An important feature of the simulation is 
that the line response is calculated by simple multiplication and 
the convolution represented by (13) need only be calculated 
once for each line and then used in each numerical computation. 
Because of this method of simulation it is possible to calculate 
efficiently the characteristics of large networks. 

Although the computation time depends upon computer speed 
and efficiency of program code, an estimate of the computation 
time may be obtained by considering the number of operations 
required to perform the convolution delay. The computation 
time required for the static boundary conditions and dynamic 
elements must be added to the convolution computation time; 
however, for systems employing primarily static boundary con- 
ditions the convolution computation time tends to be most 
significant. The total computation time required by the simula- 
tion is: 


6. a 
Te = Ne 2) m Wnt ) ° Tx0; 


t=1 i=] 


convolution computation time 

number of time steps 

number of line segments 

number of convolution steps in line segment 7 
time to compute one cycle 


Tse; time to evaluate boundaries of line 7 

For the simulations considered in this paper, the time 7',, has 
been the dominant term in (33) and is equal to approximately 
twice the time required to perform a floating point multiplication 
and addition. For minicomputers, 7; is typically 50 to 1000 
microseconds, while for large computers, 7’,, is typically 5 to 50 
microseconds. 

To determine computation time requirements, a commercial- 
ly available minicomputer* with a 0.5 microsecond cycle time 
was programmed in FORTRAN to solve general pneumatic 
transmission networks. The fluid network shown in Fig. 8 was 


*Interdata Model 80. 
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simulated with the results shown in Fig. 9. For this case, equa- 
tion (33) indicates a computation time of: 


Te = 250 (948 T,, + 6 Tso] (34) 


The program was monitored during simulations. The total 
computation time T. was 10.9 seconds, of which 9.9 seconds was 
spent in the convolved delay calculations. This time corresponds 
to a 7’, of about 46 microseconds. In another benchmark test 
it has been found that an IBM 370/165 is about 12 times faster 
than an Interdata Model 80. On this iarger computer the program 
could require about 1 second of computation time. 


Comparison of Simulation Results With Experi- 
mental Results 


To illustrate and evaluate the digital simulation method, the 
responses to a transient pulse input of the systems described herein 
have been measured experimentally and compared with com- 
puted responses. The first case considered consists of an acoustic 
driver sending a pressure pulse into a 2.44 m long by 0.436 cm 
diameter polyflow line which is blocked at the downstream end. 
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Fig. 7 Flow chart of simulation procedure 
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Fig.9 Results ot simulation of pneumatic network 
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Fig.10 Experimental and simulated pulse response of a blocked line 


The pressure pulse delivered into the line and the pressure at the 
blocked end were measured with microphones and are displayed 
in Fig. 10. 

The experimental data were recorded by an on-line computer 
at intervals of 1.1 X 10-4 seconds and stored on punched 
cards. The experimental measured input was used as the input 
to a line in the computer simulation. The experimental and 
simulation outputs and inputs were then plotted in the same 
frame for comparison. 

The experimental system was simulated using the wave 
scattering method described before with the boundary condition 
at end 1 provided by the pressure input and at end 2 by the 
blocked condition as described. In the simulation the sample 
time At = 1.1X10-‘ s and the convolution was carried out for 
500 time steps with m = 200, n = 58, and N = 43. It is noted 
that because the input pulse waveform effectively has a lowest fre- 
quency component greater than 50 Hz, a value of m appropriate 
for this range has been used. The pressure input p:, was deter- 
mined from the experimentally measured driver input pulse. 
All other parameters in the simulation were established from 
line georaetry and the properties of air at standard conditions. 
The values of N and 7’, for the line are: 


N = 43 Ta = 0.007 s 


The simulation was performed on an IBM 1130 computer and 
required 15 seconds in computation time. The results of the 
simulation the experimental measured 
As indicated in Fig. 10, the 
simulation reproduces the response with very good detail in- 


are compared with 
blocked end pressure in Fig. 10. 


cluding the delay, the initial peak pressure amplitudes, and the 
subsequent reflected pressure amplitudes and waveform. 

In a second experiment, the blocked end was replaced with an 
orifice which had a measured‘ CyAo = 0.0049 cm? and the pres- 
sure at the line input and at the orifice end of the line measured. 
A digital simulation of the experiment was also conducted using 
the same method and parameter values as described for the 
blocked line except the boundary condition at end 2 has been 
replaced by the boundary conditions for a quadratic orifice. The 
results of both the experiment and the simulation are presented 
in Fig. 11. The agreement between the simulation and the experi- 
mental data is good. The effect of the orifice is indicated by the 
lack of reflected pulses, as may be seen by comparing Fig. 10 and 
Fig. 11. 

The final system considered consisted of three fluid lines all 
of equal cross section area joined at a junction as shown in Fig. 
12. The junction was a commercially available polyflow tee 
which had been drilled out to have 0.437 cm internal diameter 


‘Determined from steady-state pressure-flow measurements. 
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Fig. 11 Experimentation and simulated response of an orifice ter- 
minated line 
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Output Termination 


Blocked 
Open 


Measured —> Orifice 


Input Acoustic 
Driver 
Length N To 

0.52m 165 0.0018sec 

118m 83 0.0041sec 

102m 95 0.0035sec 


Fig. 12 Apparatus for junction experiment 


” 


“clean through’’ in all legs. The acoustic driver was used to 
provide a pressure input to line a. Line b was blocked with the 
pressure measured at the blocked end. Line c was blocked for 
the first test, left open for the second test and terminated with 
an orifice for a third test. The values of Reynolds numbers and 
delay times for the lines are included in Fig. 12. 

A simulation of the system was conducted with each line in 
the system represented as described for the previous single line 
tests. In all the convolution integrals At = 1.2 X 1074s and m 
= 100. The results of the simulation and the experiments in 
which the input pressure to line a and the pressure at the blocked 
end of line b were measured are displayed in Figs. 13, 14, and 15 
for the three cases of line c blocked, open, and terminated with 
the same orifice as described before. The experimental and simu- 
lated data including the delay times, the initial peak pressure 
amplitudes, and the subsequent reflections from the junction and 
other lines agree very well for all three cases tested. 


Summary and Conclusions 


This paper has presented the development of an efficient time 
domain simulation method for systems consisting of linear one- 
and nonlinear, 
The method incorporates an ap- 
proximate evaluation of the transform of the line frequency 


dimensional, uniform transmission elements 


lumped, dynamic elements. 


domain representation to determine the time response and is 
designed to be efficient in terms of both computation time and 
storage requirements so that minicomputers may be used fer 
simulation. The simulation of networks consisting of line elements 
with nonlinear terminations has been performed and compared 
with experimental data. The agreement between the experi- 
mental data and the simulation is excellent. 
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Fig.14 Experimental and simulated junction structure with one open 
end 
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Planar Thermic Elements for Thermal 
s sucxtey @ Control Systems 


Associate Professor. ay ° —s 
Thermics is a newly emerging control discipline where heat flows are modulated by 


RL ART temperatures. Heat is the only energy medium necessary for control, analogous to other 
. L. Meb HY single-medium control disciplines such as electronics, hydraulics, pneumatics or 
NSF Fellow. fuidics. Although initially dominated by elements based on vaporization heat tranfser 
Department of Mecheniet Eanlacering. of the heat pipe, a new class of planar ther mic logic and control elements has been 
Massachusetts Institute of Technology, invented that relies on conduction and convection heat transfer. Three different tempera- 
Cambridge. Mass. ture modulated thermal resistors.are presented to explore the various types of modula- 
tion possible. The convection type is shown in several variations (amplifier, NOR 
gate, and diode) to illustrate analog, digital, and nonlinear control elements. Bond 
Graph models and experimental data are presented on characteristics of these proto- 

type devices. 


Introduction elements which can be combined together to give more complex 
control systems is essential to the formation of a control dis- 
cipline. 

These control elements can be loosely categorized as either 
active or passive. Common usage defines active elements as 
those which require power supplies (multiports) and is further 
broken down into analog and digital elements. On the other 


A very important category of control is thermal controi. 
Most common thermal control systems use several energy media. 
For instance, in a house-heating system an electric energy medium 
senses and controls the flow of heat energy to the house. A car’s 
coolant control system also uses two media: hydraulic energy 
provided by the water pump and heat energy in the coolant flow. 


he . ' hand passive elements are either linear (resistors, capacitors, or 
Most other control disciplines—electronics, hydraulics, pneu- 


: yee ; : inductors) or nonlinear (e.g., diodes). This classification of 
— om Suidico—cperate in wad a single energy medium. elements is admittedly vague with imprecise boundaries between 
Electronic control systems use electric power supplies and manip- 
ulate electric flows; hydraulic control systems use hydraulic 
power supplies and manipulate fluid flows. The authors examine 
the thermal area to create a control discipline based on heat 
energy. Thermics—for thermal logic—uses heat energy sources 
and modulates heat flows. 


categories. However, at least the categories are familiar to most 
control engineers. 

Typical element classification for the most common single- 
medium control discipline is shown in Table 1. Also isted are 
the associated primary variables for each discipline. Primary 


eh ‘ : variables—composed of the effort variable and the flow variable— 
Single-medium thermal control systems have existed for years. 


Indeed, one of the first evidences of feedback control was an 
alchemist’s incubator in the seventeenth century—it was in 
essence a single-medium thermal control system [1].! Other 
single-medium systems (primarily temperature controllers) are 
currently being researched [2, 3]. The basis of this discussion, 
however, is to present single-medium thermal control systems 
from an elemental viewpoint. Other control disciplines rely on 
interconnections of a few basic control elements for both signal 
flow and power flow manipulation. The development of control 


define the nature of passive elements. A linear relation between 
primary variables gives a linear element, a nonlinear relationship 
gives a nonlinear element [4]. For example, electronics has pri- 
mary variables of current and voltage, the relation between 
which define the familiar linear passive elements: resistors, 
capacitors and inductors. They also define the passive nonlinear 
elements (hereafter called simply nonlinear), for example, the 
diode. Typical electronic active elements are digital gates and 
analog amplifiers. 

The last column of Table 1 is headed thermics. The appro- 
ars priate thermic primary variables are temperature (the effort 
‘Numbers in brackets designate References at end of paper. variable) and heat flow (the flow variable). Most es 
: familiar with thermic resistors (insulators) and thermic capacitors 
Contributed by the Automatic Control Division for publication in the (masses). How sa il the 1 lecade tl . aim the - 
Jou RNAL OF Dynamic Systems, MEASUREMENT, AND ConTROL. Manuscript copes és: owever, until the last decade there were no - norman 
received at ASME Headquarters, November 12, 1976. active or nonlinear control elements. To be sure, certain non- 
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Table 1 


Hydraulics Electronics 


Primary variables 
Flow 


Effort 


Current 
Voltage 


Liquid flow 
Pressure 


Typical passive elements 
Linear (R, L, C) 
Nonlinear 


Inductor 
Diode 


Accumulator 
Check valve 


Typical active elements 


Analog Spool valve Transistor 


Digital Pressure switch Flip-flop 


(2) Discovered in last decade 


linear effects have been uncovered: glass admits more radiant 
heat from a high temperature source than from a low tempera- 
ture source; a horizontal water layer will allow more heat transfer 
from bottom to top than from top to bottom. However, these 
effects are difficult to embody into systems the way a hydraulic 
check valve or an electronic diode can be “wired’’ into a system. 


Thermic Elements 


Just as modern electronics is based on transistors and diodes, 
thermies is based on thermic amplifiers, gates and diodes. These 
thermic elements must, by definition, use only heat flow and 
temperature for sensing and modulation as well as for “power 
supplies.’’ They should also be simple elements, preferably with 
no moving parts. The fluidic control discipline emerged from 
pneumatics primarily because fluidic amplifiers have no moving 
parts. 
preferred over sliding parts; jet pipes have replaced spool valves 


However, if a moving part is required, flexing parts are 


as the first stage of many hydraulic servo-valves precisely for this 
reason. Lastly, these elements should be easily interconnected 
to form more complex control systems. 

Fig. 1 shows a field effect transistor (FET) supplied by a 
voltage source. The effort variable (voltage) at one port of the 


























Fig.1 Electrical-thermal analogy 
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Pneumatics 


Air flow 
Pressure 


Capillary 
Check valve 


Flapper nozzle 


Diaphragm valve 


Fluidics 


Fluid flow 
Pressure 


Tank 
Vortex diode 


Stream interaction 
amplifier 
Turbulence 
amplifier 


Thermics 


Heat flow 
Temperature 


Insulator 
Thermic diode® 


Thermic 
amplifier 
Thermic 
switch® 


transistor modulates the flow variable (current) at the other 
port. In this configuration the transistor acts like a voltage- 
controlled resistor. Fig. 1 also shows a conceptualization of the 
analogous active thermic resistor supplied by a temperature 
source. The effort variable (temperature) at one port modulates 
the flow variable (heat flow) at the other port of the resistor. 
The ‘‘power supply” is the temperature source. Both of these 
devices are members of a device class called “effort controlled 
resistors’? (ECR’s): an effort variable on one port modulates the 
flow through one or more other ports. 


Vaporization Devices 


Temperature-controlled thermal resistances have been known 
since the early 1960’s when the first researchers in heat pipes 
noticed that noncondensing gases could be used to vary the 
thermal resistance of a heat pipe. The noncondensing gas 
“blankets”? a portion of the heat pipe’s condenser section, re- 
ducing the amount of heat transferred by the pipe [5]. A typical 
configuration of this “variable conductance heat pipe’’ is shown 
in Fig. 2; with proper design, heat flow through the device can 
be made a strong function of the sensor volume’s temperature. 

The variable conductance heat pipe meets all the requirements 
of an active thermic resistor (ECR): simple, single-medium, 


CONTROL 


TEMPERATURE 


SENSOR 
VOLUME 


MODULATED 


HEAT FLOW 


Fig. 2 Variable conductance heat pipe 
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Fig. 3 Film ECR 


elemental form. Both active (analog and digital) and nonlinear 
elements have been developed based on the vaporization processes 
of heat pipes [6]. In the last decade, many thermal control sys- 
tems based on these elements have become operational; a small 
sample of applications are cited in the references [7, 8, 9]. 

Most thermic systems incorporating heat pipes have been 
serospace applications. Heat pipes, which usually rely on surface 
tension to convey condensate to the evaporator section, work best 
in space where surface tension is a dominate force. In a gravity 
field, heat. pipes must have their evaporator located beneath 
their condensers, often an awkward situation. Furthermore, 
fabrication and development costs of heat pipes has kept their 
price high. Hence there have been few terrestrial applications of 
control systems based on heat pipes, although a few experimental 
systems show promise [10]. 


Planar Devices 


A new class of thermic elements is being developed by the 
authors at M.I.T. 


zation heat transfer as the heat pipe devices, these elements are 


11, 12, 13). Rather than being based on vapori- 
based on convection and conduction heat transfer. Because the 
conduction and convection modes require thin components with 
large surface areas, these new devices can be classified as “planar’’ 
thermic elements. Heat pipes, on the other hand, usually have a 
long tubular structure (although planar dissipators are often 
attached to the ends for impedance matching to a surrounding 
fluid). 

Three types of planar thermic devices are shown in Figs. 3, 4, 
and 5, in an ECR configuration, with both negative and positive 
input sensors. Each device has a planar modulator associated 
with the sensors: the modulator varies its thermal resistance in 
response to the temperatures of the sensors. In each case, the 
sensor is 8 Volume of fluid which responds to a temperature change 
by expanding. This expansion is then “piped’’ to the modulator 


where it changes the modulator’s thermal resistance. Hence, 
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= as ae 
Fig.4 Pneumatic ECR 


the thermal resistance is controlled by the effort variable (tem- 
perature) of the sensor: an effort-controlled resistor or ECR. 
Note that expansion due to temperature change in the sensor can 
be either thermal expansion or the varporization of the fluid con- 
tained. 

Each type of ECR changes its, thermal resistance in a different 
manner. In the film ECR, Fig. 3, the positive sensor’s expansion 
drives a liquid film between two conducting plates; since liquids 
conduct heat better than gasses, the thermal resistance of the 
modulator is reduced as liquid rises between the plates. Using 
air and water, a ten-to-one change in thermal resistance is typical. 
In the pneumatic ECR, Fig. 4, thin diaphragms are deformed by 
the expansion of one sensor’s fluid, changing the gap between 
the diaphragms and outer plates. Thermal resistance is pro- 
portional to the average gas thickness, at least until convection 


If the 


diaphragms touch the outer plates, variation of thermal resistance 


eddys become the predominate mode of heat transfer. 


with sensor temperature is single-valued but. more complicated 
Note that the pneumatic amplifier has moving parts (of the 
flexing variety). The convection ECR, Fig. 5, using the modulated 
thermosiphon of Buckley, is more complicated because its thermal 
resistance is only indirectly modified by 


fluid [14]. 
by a thick slab of insulation. 


expansion of sensor 
It consists of two enclosed layers of liquid separated 
The vertical-oriented layers are 
connected top and bottom by ducting to form an enclosed cir- 
culation loop. W hen one lave is heated, buoy ancy forces cause 


the other layer to be heated by natural convection. Expansion 
of the negative sensor fluid is made to block this circulation loop, 
expansion of the positive sensor fluid opens the loop. One 
method is a “bubble trap’’: negative sensor fluid expansion in- 
jects a small air bubble in the trap preventing fluid circulation, 
while expansion in the positive sensor reduces the bubble’s size. 


In all ECL’s, 


sensor reduces thermal resistance through the modulator 


these a temperature increase at the positive 
while 
an increase in temperature at the negative sensor has the op- 


posite effect. Thus they operate differentially: the temperature 
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Fig.5 Convection ECR 








difference between the positive and negative sensor changes the 
device’s resistance. Other variations, such as single-input ECR’s, 
are possible but less useful. 


Convection Devices 


Variations and connections of each type ECR will produce 
amplifiers, NOR gates and diodes. These three elements are the 
most important elements of analog, digital, and nonlinear con- 
trol, respectively. The amplifier, in combination with the passive 
elements (resistor and capacitor) form virtually all types of 
analog control systems.? The building block for digital control 
systems is the NOR gate; both combinational and sequential 
logic devices can be formed from NOR elements, e.g., AND gate, 
flip-flop [5]. A rich and varied control discipline should contain 
nonlinear elements as well as analog and digital elements; a very 
basic nonlinearity is the diode since combinations of diodes and 
amplifiers can form any nonlinearity [16]. To show how the 
various planar ECR’s can be modified into the other important 
elements, one type—the convection ECR—will be presented in 
much more detail in the discussion which follows. Experimental 
results are included for each element to give a feeling of the actual 
temperatures and heat flows involved in our current research. 


Convection Amplifier 


An amplifier can be created out of two active temperature 
modulated resistors. Fig. 6 shows a schematic with two convec- 
tion ECR’s bonded to a common junction which acts as the out- 
put terminal. The ECR’s and junction are sandwiched between 
a temperature source and a temperature sink. Since each re- 
sistor has two inputs (one to increase resistance and one to de- 
crease it) they can be connected as shown to create a differential 
amplifier. 

If some arbitrary temperature is called “ground,” so that all 


*Note that only one “energy storage element”’ is required for analog control. 
Hence the thermic inductor—which doesn't exist—is not required. 
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Fig.6 Amplifier schematic 


other temperatures are measured in reference to it, then an ex- 
pression for the gain of the amplifier can be derived. Let the 
temperature inputs on the plus and minus terminal (measured 
relative to ground) combine so that: 


Tin = Tt — T- 


Then the resistors respond to this input temperature in the fol- 
lowing manner, where ry is their maximum resistance (due to 
conduction through the insulating slab) and rz is their minimum 
resistance (due to thermosiphon convection): 

ty + rh 


Ry = — = f(Tin) 


Rr wre + f(Tia) 


Letting the temperature source be the same number of degrees 
above ground temperature as the temperature sink is below 
ground, Kirchoff’s law applied to the junction between the two 
ECR’s will yield the following expression for output temperature 
in terms of input: 


T, - 7 
, h ¢ " 
Tout = ———~ f( Tia) 
ty + TL 

Where 7’, and 7’, are the source and sink temperatures, respec- 
tively. To a first order f(7'in) can be approximated by f( Tin) 
= KT in for small temperature ranges. Then the gain can be 
expressed as: 
: nak 

a K I & 
tv + TL 


GAIN 


Tout -_ 


The amplification occurs because the temperature input re- 
quired to change the resistance of 2; and Ry is small compared to 
the temperature change at the output. For a small positive tem- 
perature input, Ry becomes conductive while Rz, becomes an 
insulator. The result is the junction has good thermal com- 
munication with the heat source and poor communication with 
the sink, so its temperature moves closer to the source tempera- 
ture. The reverse happens for a negative input. The gain is the 
ratio of junction temperature change to the amount of input 
required to effect this change. 

Feedback techniques applied to amplifiers such as this lead 
to all the usual analog control systems: 
and so forth. 


controllers, lead-lags 
Development is underway of an operational 
amplifier to investigate feedback techniques of planar thermic 
amplifiers. Note that -the time constant of these devices is very 
large—on the order of minutes—-compared to electronic, hy- 
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Fig. 7 Bond graph of ECR 


INSULATION 


draulic, or fluidic control systems. However, thermics is in- 
herently suited to the medium which it is controlling: thermal 
phenomena usually occur slowly because of the high thermal 
capacitance associated with most substances and the high thermal 
resistance associated with convection heat transfer to air. 

To predict the dynamic thermal behavior of the amplifier it 
can be considered a dynamic system and modeled appropriately . 
The building block of the ECR is the convection resistor, which 
ean be characterized by the bond graph shown in Fig 7, using 
temperature/heat flow as the effort/flow variables. This basic 
bond graph structure performs many different functions de- 
pending on how the flow resistance, Rr, is modulated. For the 
case of an ECR the resistance is modulated by the effort variable 
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Fig.8 Experimental convection amplifier 
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of a separate sensor; two such ECR’s are combined to give an 
amplifier. The flow resistance, Rr, can be modulated other ways 
to give other devices that will be discussed later. 


Convection ECR 


In the experimental convection ECR, the sensors consisted of 
thin metal bags containing a small amount of Freon 1ll—a 
fluorocarbon liquid whose boiling temperature is about 24°C 
(75°F) at atmospheric pressure. The bags were enclosed by flat 
polycarbonate boxes whose only opening was connected by vinyl 
tubing to the bubble trap at the top of the heat flow modulator. 
The modulator itself consisted of two open-topped metal con- 
tainers separated from each other by a Plexiglass frame. A 
block of foam inside the frame formed a liquid layer next to each 
container; a passage at the top and a bubble trap at the bottom 
completed the circulation loop (see Fig. 8). 

For testing purposes, one container was filled with boiling 
water while the other was filled with an ice slurry. Hence the 
temperature drop across the resistor was nearly a constant 109°C 
(180°F). The sensor was suspended in a temperature-controlled 
liquid bath to vary its temperature; heat flux gages (measuring 
temperature drop across a fixed thermal resistance) determined 
the heat flux through the amplifier. During operation, an in- 
crease in temperature of the sensor caused vaporization of the 
Freon, forcing air into the modulator’s bubble trap. Stability is 
aided by pressure feedback: as the Freon vaporizes, the slight 
spring rate of the metal bag increases the vapor pressure on the 
remaining liquid and prevents further vaporization. Fig. 9 
shows the experimental results by one of the authors (McCarthy) 
of heat flow through the modulator versus temperature of the 
negative sensor. 

Since the heat flow measurements were taken across a constant 
temperature drop, the graph indicates the variation in thermal 
resistance; total change in thermal resistance is about 12 to 1. 
As shown previously, an amplifier is built from two such re- 
sistors connected through a conducting junction block. Because 
a temperature input of only 1.5°C (3°F) will drive the resistors 
from saturation to saturation, an amplifier constructed from 
two of them would have a single stage gain of 48. 
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Fig. 9 Convection amplifier characteristics 
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Fig. 10 Thermic NOR gate 


Convection NOR Gate 


Since the convection ECR saturates, it can be used in a digital 
mode. A NOR gate operates in the following fashion: it is always 
ON except when any input is ON (in which case the gate turns 
OFF). In the experimental device, ON was chosen as a tem- 
perature higher than ambient and OFF was chosen as ambient 
temperature itself. Hence the output surface of the gate will be 
hot unless any input is hot. 

The experimental device consisted of a ECR 
modified to have two “bubble traps’’ in series, each with its own 


convection 


sensor. Referring again to Fig. 7, the flow resistance, Rr, of the 
bond graph strueture was modulated in a digital manner by the 
effort variables of both sensors. A constant heat flow was applied 
to one surface of the gate’s modulator by radiation heat transfer. 
The temperature of the other plate (the gate’s output) was 
monitored with thermocouples. Temperature changes were im- 
posed on sensors by emersing them in a beaker of water of a 
particular temperature. Fig. 10 shows the average temperature 
of the output plate for a typical input combination. With both 
sensors at ambient, the temperature of the output plate was much 
hotter than ambient. However, with either or both sensors hot, 
the output plate cooled nearly to ambient. Like the convection 
amplifier, the time constant of the NOR gate was about 15 
minutes. Further development work is in progress by the authors 
investigating interconnections of NOR elements; specifically, 
fan-in and fan-out 
flip-flop. 


characteristics and the construction of a 


Convection Diode 


The convection ECR can also be modified into a diode. A 
thermic diode transmits heat easily in one direction but poorly 
in the opposite direction; its thermal resistance differs depending 
on which direction heat is being transmitted. The modification 
involves placing a check valve at some point in the circulation 
loop so the natural convection flow will be blocked in one diree- 
tion but will be free to flow in the other direction. 


The check valve must be a very special one: it must respond 
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to very slight pressures, it should seal completely in the back 
flow direction, and it should preferably have no moving parts. 
Fig. 11 shows the check valve of Buckley which meets these 
requirements. A riser tube juts slightly through the interface 
of an oil layer floating on water. If the pressure in the riser is 
higher than the liquid surrounding it, water will flow up the 
tube and out over the tube’s end. In experimental models de- 
veloped by Buckley only 6.9 Pa (0.00L psi) pressure was re- 
quired to initiate flow; the check valve owes its sensitivity tothe 
buoyancy-caneelling effect of a two-liquid manometer. When the 
reverse pressure gradient is applied to the check valve, oil is 
drawn down into the riser tube until a pressure balance is 
achieved—the buoyancy difference between oil and water cancel 
the applied pressure—and flow ceases. The overall effect is a 
check valve which responds to very slight forward-flow pressures 
yet seals completely in the back-flow direction. 

This check valve installed in the circulation loop of a convec- 
tion modulator produces a thermic diode. Forward heat flow is 
determined by convection heat transfer of the circulation loop; 
reverse direction heat flow is related to the insulation thickness 
between the liquid layers. Referring to the bond graph of Fig. 7, 
the flow resistance, Rr, is made a nonlinear function of the tem- 
perature difference between the two sides of the resistor. In the 
experimental thermic diode, one side was heated electrically 
while the average surface temperature of both sides was measured. 
Fig. 12 shows the heat flow through the diode as a function of 
temperature drop for both forward and reverse bias. For this 
configuration, having a planar area of about 0.56 m? (6 ft?), the 
diode had a thermal resistance in the reverse direction about 40 
times that in the forward direction. No-flow heat transfer was 
dominated by conduction heat flow through 0.076 m (3 in.) of 
closed cell foam installation. 

Thermic diodes are being developed by the senior author for 
solar energy applications [17]. Diodes are constructed in 1.22 m 
X 2.44 m (4 ft X 8 ft) panels which replace the roof of a building. 
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Fig. 12 Thermic diode characteristics 


When the sun heats up the outer surface, solar heat. is trans- 
mitted inside, at night or on cloudy days the diode prevents 
heat loss from the building. In one variation, the inner water 


layer is enlarged to provide a heat storage capability as well. 


Summary 


Thermics is a single medium control discipline based on simple 
control elements such as amplifiers, gates and diodes. Like other 
control disciplines, these elements can be interconnected to form 
complex control systems. Thermies uses temperature to control 
heat flows; heat sources provide the only energy to thermic sys- 
tems. 

Early research in thermics evolved around vaporization de- 
vices, those based on heat-pipe principles. More recently, a new 
class of thermic elements shows promise for terrestrial applica- 
tions of thermal control. Three types of these devices—called 
planar thermic elements—are being developed by the authors: 


film devices, pneumatic devices and convection devices. Planar 
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elements are distinguished by a panel-like geometry which pro- 
vides efficient heat transfer characteristics. 

Each type of device can be modified into the basic thermic 
elements (amplifier, NOR gate and diode). For illustration, these 
variations are shown for the convection device along with relevant 
experimental data. 
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Professor, The motion of an eccentric rotating shaft and disk is studied for an exponential transi- 
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University of Notre Dame, is found for the response, and it is shown that this response can yield higher amplitudes 


Notre Dame, Ind. when compared to previous work where the angular velocity is varied linearly. 


Introduction 


The response of a rotating shaft and disk with an eccentricity 
was analyzed in Reference [1]? and summarized in Reference [2]. 
In these works a rigid disk is attached to a weightless elastic 
shaft supported by rigid bearings as shown in Fig. 1. The mass 
center of the disk is assumed to be eccentric from the geometric Define 
center by the distance, a, as shown in Fig. 2. These analyses 
assumed that the shaft’s angular velocity changed (increased or 
decreased) linearly with time and passed through the critical 
speed of the shaft-disk system. Generally it was found that the 
peak displacement response while passing through resonance was 
less than the corresponding steady state resonant amplitude. 
Furthermore, this peak response lagged the time at which the 
shaft’s speed equaled the critical speed. 

When a shaft is sped up through a critical speed by a constant- 
torque motor or is freely slowed down, the speed change is more 
realistically exponential than linear. This is particularly true 
when the final rotating speed is close to the critical speed. It is 
shown in the following analysis that the response can be found 
in closed form when the angular velocity varies exponentially. 
This analysis is based on a solution obtained for a similar prob- 
lem dealing with the passage through resonance of a free flight 
missile [3]. 





General Solution 


The equations of motion for the system in Figs. 1 and 2 are 
given in Reference [1], for the case with linear damping: 


'This work is based upon a portion of the PhD dissertation submitted by 
B. M. Naveh to the University of Notre Dame. 


2Numbers in brackets designate References at end of paper. 
Contributed by the Automatic Control Division for publication in the 


JournaL oF Dynamic SysTeMs, MEASUREMENT AND CONTROL. Manuscript 
received at ASME Headquarters, September 16, 1976. Fig.1 Shaft and disk 


48 / MARCH 1977 Transactions of the ASME 


























Y* 


Fig. 2 Disk geometry and coordinates 


c= r+ ty (3) 
Multiply equation (2) by 7 and add to equation (1) to get: 
ee ka 
#4+— 2+ z= e'? (4) 
m m m 
The prescribed exponential angular velocity change is given by: 
d = dee + (1 — eM) (5) 


where @., is the final steady state angular velocity and go is the 
initial angular velocity. 

For the increasing angular velocity case (i.e., 6 > 0, d. > wo 
> do and L < 0) the solution is given by 


idr ipr 
®,{ 1, a + 1; - ®,{ 1, pe t+ 1; 
wat ( . L ) ( L ) e*? 


2 - ibe + Te 
+ Kyes' + Kees" (6) 


z(t) = 


fi = bn — 


i(do * wo) — A 
M,2: >= I 


nas = 


kn - 


m 


Nomenclature 


wo = y: 
m 
®,, are the degenerate hypergeometric functions [4]: 
<a co L 
t 
hs (nmett B)a14+ De a 


"TIT (un; 2 + 9) 


q=l 
K,,2 are determined by the initial conditions, 


_ t—-tm. G heme — ido) — %,(0) 
K,; = ——— + — 





2iwy a _ ide + Ti 


$,(0)(r2 — ido) — (0) 
ee) on 





. 1120 — Zo 
K,; = —_— - 


a [2:0)(r1 — ide) — ,(0) 
2iwo 4 


ido tn 





$,(0)(r1 — ido) = #,(0) (14) 
= ideo + T: 

The ®,, are converging series and can easily be computed 
on a digital computer. When the shaft decelerates through the 
critical speed, @ < 0, L < O and do > wo > dw. The correspond- 
ing general solution is obtained in the same manner as above. 
This solution is given in detail in Reference [3]. There it is also 
shown that gyroscopic coupling can also be included in equation 
(4), if desired. 

As done in Reference [1], it is convenient to define a nondi- 
mensional amplification factor, RM, for the response. This is 


RM(t) = 


z(t) 
a 


(15) 


The peak value of RM as the system passes through the critical 
speed (linearly with time) was found in [1] to be lower than 
when the system is in a sustained resonant condition. When 
the passage is exponential the results can differ significantly as 
indicated by the following example. 


Computation and Results 


From Reference [3] the amplification factor decreases with the 
increasing of @. for a given ¢o. The linear law for the angular 
velocity can be regarded as a limiting case of the exponential law 
where ¢, — ©. Therefore the amplification factors obtained in 
Reference [1] are low for a design where @¢,. is close to wo. To 
demonstrate this, the amplification factor for both cases will be 
compared using the example given in reference [1]. The ex- 
ponential transition is made to yield the same angular accelera- 





= eccentricity of disk 
linear damping coefficient 
v-1 
constants 
= transverse bending stiffness of 
shaft 
= angular velocity damping factor Mi, 
mass 
n,q = integers 
RM = amplification factor 


time 
=2+ wy 
—c/m 


spin angle 


Po — P 
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+ iwo + Xr 
= transverse, orthogonal displace- (") 
ments of the disk “7 


lide F wo) — AI/L 


= degenerate hypergeometric func- © 


tion 
Wo = (k/m)! 
derivation with respect to time 
| | = absolute value 
(0) = att =0 


Subscripts 


R = at critical speed (resonance) 
0 at t 
att 
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Fig. 3 The response of a rotating shaft no friction 





1.20 


tion, dr, through the critical speed, wo, as the example and ¢~ 
is taken to be 115 rad/sec. For. these values, the decay constant 
for the angular motion is: 


or 


L= ; - 
wo — Do f 5 


= — | sec"! 


(16) 


The results are compared in Fig. 3, for the frictionless case and 
in Fig. 4 for the friction case, where c > 0. 

In Figs. 3 and 4 it can be noticed that the maximum response 
lags the occurrence of @/wo = 1. In the case of the exponential 
variation of the angular velocity the maximum response will 
generally lag less then in the case of the linear variation [3]. 

In the linear variation case the amplification factor is the same 
whether the shaft accelerates or decelerates through resonance. 
It is shown in Reference [3] that for the exponential variation of 
the angular velocity the amplification factcr when speeding up 
is not generally the same one that is obtained when the shaft is 
slowed down, although |r| in both cases are made to be equal. 
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Conclusions 


The response of a shaft accelerating through a critical speed is 
given for an exponential angular velocity variation. Previous 
results were for a linear variation. It was found that the linear 
change of angular velocity yields optimistic results in comparison 
to those obtained for an exponential change of the angular veloci- 
ty. 


References 


1 Dimentberg, F. M., Flerural Vibrations of Rotating Shafts, 
Production Engineering Research Association, Butterworth 
London, 1961, pp. 42-90. j 

2 Loewy, R. G., and Piarulli, V. J., Dynamics of Rotating 
Shafts, Shock and Vibration Information Center, Code 8404 
Naval Research Laboratory, Washington, D.C., 1969, pp. 45-52. 

3 Naveh, Ben Zion M., “On the Dynamic Response of Mis- 
siles with Varying Roll Rates,’ University of Notre Dame, PhD 
dissertation, 1975. 

4 Gradshteyn, I. S., and Ryzhik, Tables of Integrals Series and 
Products, Academic Press, New York and London, 1965, Fourth 
Edition, p. 1058. 


Transactions of the ASME 





State Estimation and Parameter Identification 
of Freely Spinning Flexible Spacecraft 


This paper studies the attitude estimation of a flexible spacecraft using the noise-cor- 
rupted measurements provided by sensors. The flexibility of the system is expressed in 
terms of internal deformation modes, and only a limited number of the latter are kept 
in view of simplifying the model. The estimation method is based on the Extended 
Kalman Filter algorithm, using an augmented state vector which includes the unknown 
parameters of the system. The filter thus provides an estimation of these parameters 
along with that of the original state-variables of the system. The method is implement- 
ed using a simulated model, thereby permitting an appraisal of the results, and refer- 
ence is made to results obtained for a practical case where a somewhat different approach 
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was taken as regards the way in which the flexibility is expressed. 


Introduction 


The estimation and identification problem has already been 
studied for rigid spacecraft in a previous paper [1].!_ The intro- 
duction of nonrigidity however substantially alters the approach 
to the problem [2]. If the same model were kept, most param- 
eters would become time-varying and it was shown [1] that this 
is prohibitive to the proper functioning of the filtering algorithm. 
This led us to change the model of a general system by intro- 
ducing all deformations through two basic low-frequency modes 
of the non-rotating system, the parameters to be identified then 
being the coefficients of the amplitudes of these modes. After 
deriving the equations of motion based on this assumption (Sec- 
tion 1), the state equations of the system may be written (Section 
2). The basic theory of the filtering algorithm (nonlinear Ex- 
tended Kalman Filter) is ‘hen outlined in Section 3, and the fol- 
lowing sections give details concerning the application of the 
algorithm to the problem at hand. The numerical results of an 
application to a simulated model of SKYLAB are given in the 
penultimate section as well as some guidelines as to the optimal 
choice of initial values. Recent results obtained by applying 
the filtering algorithm to a model of European Satellite GEOS 
IV are reported in the last section; for this specific case in which 
the system is particularly well known and the dynamics have 
been thoroughly analyzed, real physical parameters are chosen 
rather than amplitudes of the internal deformation modes—this 
somewhat simplifies the modelling and leads to excellent results. 


'Numbers in brackets designate References at end of paper. 
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1 Dynamical Equations of the System 


The Euler-Liouville-Resal equations for a field-free spacecraft 
read: 


aXdcet+seut+5-8+h+uXh=0 (1) 


where J = inertia tensor 


h = 


internal angular momentum 


o = local derivative. 


We shall express this equation in a mean reference frame {X.} 
(extended Buckens’ frame [2]) defined in such a way that the 
relative angular momentum is equal to zero in linear approxima- 
tion; @ is then the angular velocity of this frame relative to 
inertial space: 


“ Oi 
@ = w"[Xal, We 
Wa 


a = Il, 2, 3. 


In linear approximation, equation (1) then becomes (in matricial 
form): 


aw + Jw+ Ja = 0, 


where 


— Ws Wi 


Equation (2) will yield three scalar equations for the variables 
Gi, We, Ws. 
Flexibility will be considered from a modal viewpoint. 


A de- 


formation u may be written as a sum of modes: 
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u(t) = Y).Biltos, (3) 


i-1 


where @, is the eigenvector and £;(t) the amplitude of the 7th 
mode. The n normal modes of the nonrotating system are chosen 
to form the basic set of deformation modes. The n amplitudes 
(n is the order of the nonrotating system) are the components of 
a vector B: 


B = [B, ..-, Bal’. 


The deformation equations result from classical Lagrange equa- 
tions in which the kinetic energy may be written: 

: 1 2 : 

T = F ww + wh + 3 BTM .£B, (4) 


and the potential energy of deformation will be supposed of the 
form 


1 
Ua = 3 BK.B + FB, (5) 


where F corresponds to generalized prestresses and Ka and Ma 
are (n X n) matrices. 


The dissipation function will be written 
1 ° ° 
R= - BT SB. 


In the reference frame we have chosen, 4 is diagonal and 
will be set equal to the unit matrix EF by normalizing each 
eigenvector @; in the following manner: 


di*M adi = 1, 


(d;* is the transposed conjugate of ¢;). 

This definition amounts to multiplying 6;(t) by a certain con- 
stant factor and will greatly simplify matters further on. 

As we shall be considering only small perturbations around an 
equilibrium, we will retain only linear forms of the variables w 
and @ in the equations. 

At equilibrium: 


aa 


Wm=m = 8 =0, t=1,... 
Ws = Wao. 


In order to obtain exact linearized equations, the following 
forms will be adopted for the elements of J and h in order to 
have a quadratic expression for 7’ in B;, wi, we, ws', (Ws = Wo 
+ w,'): 


Ju = 
Jn 
Ji = 0, 
Ji3 = A,7£, 
Ja = A,78, 
Jes = Ip + Az™B + BTx'B where ' = z'7, 
hs = BTCB. 
The linearized equations are then, in the frame we have chosen: 
Tye, — (Ip — Ts)ws@2 + waoAi7B — ws®As™B = 0, 
Tang — (Is — Ti)eos0e + croha?B + wtAr™B = 0, 
T 3G, + wsoAs7B = 0 whichyields w,;' = — onl? (6) 
3 
EB + (S + wal")B + MaB — wrodionr — waoArwe2 — wirAseos' = 0, 


where Ile = Ka — IIL'wye?, 
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F = Asso? 


r =ct-cC=-I™. 


(equilibrium condition), 


The last equation may be rewritten as follows, after introducing 
the expression of w,' resulting from the third equation: 


T 





: 2 \ 
EB + (S + wal)B + (1. + - Asset) — woAliar 
3 


- wsoA2we = 0. (7) 


It would be possible to rewrite the Euler equations in terms of 
attitude angles 6;: this doubles the number of rotational equa- 
tions but is necessary if the 6;’s are the variables of interest; it 
will not be done here. 

We now make the assumption that there are at most two 
dominant deformation modes in the system. This means that a 
deformation u may be sufficiently well approximatec by a com- 
bination of two modes: 


w = Bilt)or + Bo(t)pr. 


This hypothesis is quite realistic as there will rarely be more 
than two high-amplitude low-frequency modes, and al! high- 
frequency modes can exist only for a short time due to the damp- 
ing in the system. The truncation of the series (3) may result 
in a reorientation of the two remaining modes and a loss of 
significance of these modes, as well as of parameters A, II, S, 
ete. 
equations describe the dynamics of the system adequately and 
enable us to estimate the kinematic variables w;; the identifica- 
tion of the parameters is only a step towards a better estimation 
of these variables. 


This is, however, not of great importance as long as the 


We now rewrite the two deformation equations: 


By + iB, + wroloBe + 16: + ITB. —_ woAna = WsoA 1262 = 
Be + C232 + wrol'o; + Inf; + T1282 = 


wWoAnw: — warAn2wW: = 


where 


re Ee 
| -Ts oy; 
PAs 

me Aga ‘ 


(i = 1, 2). 





It should be noted that all the coefficients of these two equations 
have had their physical dimensions divided by that of Ma (in- 
ertia). 


2 State-Vector and System Equations 


Each of the parameters appearing in the dynamical equations 
may be decomposed into a nominal value and a bias, e.g., 


An : 
Au? ] 


The ratios of the bias to the nominal value will become a state 
variable y; of the system (y; = 0 if the corresponding parameter 
is nonexistent). There are, of course, alternate ways of writing 
the relation between a parameter and the corresponding state 


/ 
An = An® + An = An 4 + 
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variable, e.g., yi = parameter. The decomposition above was 
chosen by analogy with that of a previous study [1] and has the 
advantage of equating the nominal value of any parameter to 
zero. 

Moreover, the elements of the transformation matrix L, de- 
fined further on will also have to be identified. In our case this 
means 4 extra variables. 

We thus have the following state-vector (dimensionless com- 
ponents): 


"1 ’ ! ! Bi 


Wao 
Ay® A 2° ; 
tu 
TT,,°’ 
x I 
Ty? ZL * 


TT,° ’ 


ts la log 


La’ “ale” ie” 


Yie ’ 
Ln° 


The equations of the system may be written as a set 18 first- 


order nonlinear differential equations: 
y = fly). (9) 


The developed form of these equations may be found in {3}. 
They are linearized, around the equilibrium state, in variables 
yi through y.. 


3 Estimation and Identification 


The object of the study is the estimation of the basic state 
variables @, @2, 81, 82, B:, B2 describing the kinematics of the 
satellite relative to its reference frame. The data available are 
provided by m measurements taken at constant intervals in 
time, and which we assume to be corrupted by white noise: 


z= 
2 


h(y) + », (10) 


where E[vv7] = 


R (measurement noise covariance matrix) 


z= 


measurement vector 


h = nonlinear function of the augmented set of 
state variables 


v = white noise vector. 

We are therefore confronted with a classical filtering problem 
for an augmented state-vector. As the state equations are in- 
tegrated between steps, a good knowledge of the parameters ap- 
pearing in these equations is necessary, which explains why 
these parameters must be identified along with the original state 
variables of the system. 

We shall make use of a nonlinear discrete-continuous extended 
Kalman filter derived by Jazwinsky [4], and which reduces to 
the following steps. 


Between Measurements. Integration of the differential equa- 
tions describing the variations of the state estimate j and of the 
error probability covariance matrix P: 


y = fi) corresponding to y = f(y) + w, 


P = FP + PFT+Q, (11) 


where 
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P = El(y — y)\(y — 9)", (n X n) matrix, 


of 
& y 


w system noise vector, 


Q 

The elements of F are given in detail in [3]. The initial values 
of 7 and P are given at each step by the reset values of the pre- 
ceding step. 

The starting values of P and Q at the first step depend on the 
characteristics of the system. 

The result of the integration is y(—), P(—), where the (—) 
means “before reset.’’ 


(n X n) matrix, 


Elw w] , (n X n) matrix. 


Reset Equations 
w+) = 9(—) + P(—)HTYHz — ACH(—))), 


P(+) = P(—) — P(—)H*Y“HP(—), (12) 


H = ( S ) (m X n) matrix, 
ov Jy, 


Y = HP(—)H™ +R, (m X m) matrix; 


y(+) and P(+) are the initial values for the next integration. 
#(+) is the optimal estimate of y, at that time. As the iden- 
tification of the parameters progresses, #(+) will become a 
narrower estimate of y. Fort + ©, both the parameters and the 
original state variables should be correctly estimated (with er- 
rors due to v and w). 


4 Implementation and Model Simulation 


A computer programme has been written to carry out the 
various steps of the filtering procedure; the inputs are the initial 
values of the state variables and matrices P, Q, R, the nominal 
values of the parameters, the expressions of function z = h(y) 
and matrix H, and the physical data of the system. 

In a practical application, the measurements would be pro- 
vided directly by the sensors; in our case, however, we had to 
generate these measurements from a simulated model which also 
enabled us to check the quality of the estimation. 

The model chosen is a simplified model of SKYLAB (Fig. 3) 
[5]. Resolution of the dynamical equations of this system, ob- 
tained from a multi-body formalism, yields the state variables 
at any given time, in particular at the sampling moments. These 
state variables will be used to generate the measurements. 

The deformation variables represent boom deflections along 
perpendicular directions. The solar panels of SK YLAB are as- 
sumed to be rigid. 


5 Measurements and Noise Generation 


We assume that an instrument fixed to the main structure can 
measure the angular rates of this main structure about two axes 
which at equilibrium coincide with X, and X:. If the deforma- 
tions are small, the sensed rates, m, and m, are equal to 


m = @ + v1 — Ways, 
™m, = Ws + 2 + wi, 


where y; and 72 are the angles describing the orientation of the 
instrument axes with respect to the {Xq} frame. 

It is clear that for small deformations the angles y: and 7: 
are linear combinations of the state variables 8; and #:, 


Bi bers &? 


MARCH 1977 / 53 





The four elements Li; will have to be identified along with the 
other parameters of the system. ; ; 
Moreover, we assume that the deformation rates 8, and #: 
can also be measured. 
We therefore have a (4 X 1) measurement vector z 


z= hy) +», 
wrolts + Ll + yrs)ys + Ln°(1 + tr0)ye 
— Lal + yr)ys — Lal + yre)yl, 
rola + Ln + yrr)ys + Lal + yis)ys 
+ Ly%(1 + yis)ys + Ln°(1 + yr0)ysl, 
ha(y) = wroys, 


li(y) = 


hi(y) = 


he(y) = wroye. 
This yields the H matrix in a straightforward manner from 


Oh; 
H;= —. 
: Oy; 


As for the noise vector v, we shall assume it to be white noise 
with uncorrelated components, zero mean and variances given 
by the following relation: 

ae hu 
oC" = , » 
At 
where At is the numerical value of the sampling interval. This 
is an aproximate formula provided by the extension of the con- 
tinuous Kalman filter to discrete sampling. In fact Room = 
Raiscrete * At strictly for At — 0. 

The actual measurements will therefore be obtained by 
adding random variables v; (i = 1, ..., 4), generated by an ap- 
propriate subroutine, to the expressions h;(y). 

The components of the system noise vector are generated in a 
similar way, the variances being provided by 


a? = Qi At. 


6 Choice of Initial Conditions 


As only two modes are taken into consideration in the mod- 
elling of the system (low-frequency modes as explained earlier), 
it would be unadvisable to excite the higher frequency modes 
unnecessarily by an unconsiderate choice of the initial condi- 
tions. This remark is true, not only for the basic state variables 
yi... Ys, but also for the corresponding block of P, whose equa- 
tion is 


Pr, = FusP ss + PF oe™ + F upP op* + Pool on? -+- Qee, (13) 


where 
Pa Pop 
> yi > 
P., en 
Put Pont 
a Cy Fem 
pe Ppp 
Fate Pics 


Pom 
Pas 
Pe m 


Part Xs 
Pate mS 
Pom 4 x 4 


Fates 
r. Pop i: 8X8 
I 


Finm 





Mmm 2 4X 4, Fem = Pune = 0. 


We see that the dynamical equations, which may be written 
Y, = FeYs + FY, (for the first 6 equations), have the same 
eigenmodes as (13). 

This observation, together with the experience gained from 
various unsuccessful runs, led us to the following guidelines in 
the choice of the initial values: 


x(0) = (initial value of model state vector) : near to 
low frequency mode of F,,, evaluated with the 
parameters at their real value—this will gen- 


erally be the case as high frequency modes will 
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have been damped out 


corresponding low-frequency mode of F,, eval- 
uated with the parameters at their starting 
values—this mode is caliby ated by means of a 
measurement 


Xa Xa r 
= ( — ) where xo = #1 ...0(0). 


n n 


P,.(0) 


n= — : —— (from R). 
error probability 
(Ax due to Ap;) « Ap; 
Ax is obtained by a sensitivity analysis of the 
eigenmodes 
Ap; is the estimated probability of error on 
p;(jth parameter) 


Praygth column) > 


diagonal matrix with elements (Ap;)? 


diagonal matrix with elements (Am,)? 
Am, is the estimated probability of error on mx 
(kth element of L) 


Pup = 0 


= 0—except perhaps for equations (5) and (6) : 
Qs ~ 0 
Qe ~ 0 


Yz...18(0) = O as the nominal value represents the best 
knowledge of a given parameter at = 0— 
in our runs, however, some of the parameters 
were left at their real values, to reduce the 


number of parameters to be identified 


Concerning the choice of 7:....(0), it may be explained as 
follows : if one accepts the fact that high-frequency modes are 
damped out, the system should be in a state near to the low- 
frequency mode, and the initial estimate is chosen equal to the 
best estimate of this mode, i.e., with the parameters at their 
starting values; a calibration of this mode may be obtained by 
an initial set of measurements. This procedure is however quite 
stringent, and it may be disregarded for systems whose dynamics 
are well known: this is the case for European satellite GEOS IV 
(section 9) for which zero initial conditions where chosen for 
ji...5 without affecting the precision of the ensuing results. 
Moreover, the measurements used for this satellite model where 
the actual ones scheduled and not elastic bending rates which 
were chosen only as an example, albeit possibly unrealistic, in 
the general model. 


7 Numerical Results 


Various tests were first run to check the validity of the two- 
mode model as compared to the simulation model. The results 
were excellent, thereby confirming, for this specific case at least, 
the validity of our assumption. We then proceeded to identify 
one parameter using three measurements : 21, 22 and 24. 

All the parameters were set to their real value except ys which 
was given a zero initial value instead of 0.1115 (the numerical 
data of the model may be found in [3]). The evolution of ys(+) 
versus time (sampling period : 2 sec = (W30/5)) is recorded in 
Fig. 1; it may be seen that #s(+) reaches the vicinity of its real 
value very quickly (3 steps) and then stays in a small region 
around this value. The evolution of Pss speaks for itself: 


Ps3(0) = 4.107? Pas(t = 100) = 2.10-4. 

As for the 6 main variables, they are very well estimated as 
may be concluded from the following table which records the 
difference between the real and estimated values, divided by the 
maximum value, for the last few steps: 
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Fig.1 identification of one parameter 











—e ee ae . 
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Fig.2 Identification of two parameters 


: 0.44 percent measurements (this yields better results than using only 3 
measurements as one of the parameters is otherwise practically 


2 : 0.1 percent unobservable). 


ys : 0.49 percent All the parameters were set to their real values except 

y« : 2 percent ys(0) = 0 instead of 0.1115 

ys : 0.2 percent (Qs; = 10-6) tu(0) = 0 instead of —0.05777. 

ve 0.1 percent (Qe = 107%) The evolution of is(+ ) and #(+ ) versus time (sampling period : 


2 sec) is recorded in Fig. 2. The final values of js(+) and gu(+) 
The next step was to identify two parameters, using all four are 
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ys(+) ~ 0.102 
gu(+) ~ —0.059 
which is quite satisfactory. 


The table of relative errors on the 6 main state variables gives 
the following results : 


yi : 0.1 percent 
y2 : 0.7 percent 
ys : 0.6 percent 
ya : 1.5 percent 
10-*) 
10-°) 


ys :2.5 percent (Qss 
Ye: 1.9 percent (Qes 
For Qss = Qse = 0 = Qii, we obtained the following : 
s(+)tin & 0.104 (6.6 percent) 
Yu(+)tin & —0.058 (0.7 pertent) 
and relative errors on 
: 0.3 percent 
: 0.2 percent 
: 0.2 percent 
: 0.1 percent 
: 0.3 percent 
: 0.01 percent 


8 Comments on the Model and the Influence of 
Certain Factors 


The choice of SKYLAB was dictated mainly by practical con- 
siderations, as all the data were available and the model had been 
dynamically analyzed. 
this spacecraft was available. 


Moreover, a simulation program for 
For identification purposes, how- 


ever, it was not the happiest choice for two main reasons : 


1 The transformation matrix L has very small elements, 
which means that important deformations of the booms cause 
only small deviations of the mean reference frame. Moreover, 
these elements are practically unobservable and cannot, there- 
fore, be included in the parameters to be identified. 

2 The system is very close to instability due to the small 
difference between J; and J. This means that one cannot tolerate 
too great a deviation on certain parameters. We encountered 
this difficulty for ys, which, when greater than 0.06, caused an 
instability of the system. 


There remained in fact only 4 parameters which could be 
identified. 

The choice of th. initial values x(0), 7(0) and P(O) is of prime 
importance, expecially when the damping is small. In particular, 
great care must be taken in choosing P,,(0) corresponding to the 
parameters to be identified. 

As for the damping coefficients, their value is not of great im- 
portance, as long as they are sufficient to damp out the higher 
modes. This is why we did not include them in the state vector. 


‘. Results Obtained for European Satellite GEOS 


GEOS IV is a geostationary satellite (to be launched in 1977) 
consisting of an axisymmetrical rigid main body with two 
flexible cables attached to the main body at diametrically opposed 
points and carrying tip masses (Fig. 4). 
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Fig. 3 Model of SKYLAB 


Fig.4 Model of European Satellite GEOS IV 


A full dynamical! analysis has been carried out for this satel- 
lite [6] and the great confidence acquired in the model as a result 
of comparative tests led us to choosing specific physical param- 
eters instead of the more general parameters used in the modal 
analysis. 

The 13 resulting parameters include elements of inertia tensor, 
mass of main body, lengths of cables, position of mass center of 
main body, spin rate. 

The estimation/identification algorithm remains unchanged; 
measurements are provided by accelerometers and magnetom- 
eters located on the main body. 

A complete description of the results obtained may be found 
in [7]; we will limit ourselves here to a very restricted sample of 
these results which should however give a good idea of their 
quality. 


1 Identification of inertias /;, /; and Length of (Symmetrical) Cables. 
The initial (nominal) values of these parameters are 
I; = I, 
I; 
Ir-l, 
I; 


: 71.93 kgm? a= = 0.1771 


I; : 146.38 kgm? 0.5571 


T,: 20m / = 1.095 


The exact (simulated) values and the estimates obtained after 
about 30 filtering steps of 1'' each (computation time for one 
step ~ 12 sec) are as follows 


I; I; l, Qa B Y 
75.98 153.38 19.00 .2600 0.5795 1.148 
153.37 1.148 


Exact 


Estimated 76.16 18.86 .2591 0.5778 
The last. row of figures is obtained by averaging over the last 
three filtering steps, as small variations still occur (Fig. 5). 
As for the main state variables 6,, 02, and boom deflection 


angles, relative errors of about 2.5 percent are typical. 


2 Identification of Lengths of Cables. With (/:)nom. = (l2)nom- 


= 20 mand Al; = i; — (li)aom 
Al, + Al, 
hae Al, — Al 


the following results were obtained 
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Fig. 5 


las 
Exact 1 


Estimated 0.96 


3 Identification of Lengths of Cables and Spin Rate. 
(li)nom = 20 m and (wWo)nom = 1.0 
the following results were obtained: 


With 


Wo l ‘ les 


Exact “1.05 16 2 


Estimated 1.0497 15.93 1.93 


In all the tests run on GEOS, as well as on SKYLAB, the 
algorithm was still converging when the runs were brought to an 
end, as can be seen from Fig. 5. Moreover, the errors on the state 
variables and parameters were always compatible with the 
statistics, i.e., the mean deviation given by /P;; was always of 
the order (and usually greater) than |z; — 2;|. This is an im- 
portant result as it enables the user to evaluate the error by in- 
spection of Pj;. 


Conclusion 


The purpose of this study was to analyze the applicability of 
nonlinear Kalman filtering as an identification scheme for the 
parameters of a force-free deformable spacecraft, the final ob- 
jective being to estimate certain state variables of the satellite 
from a series of noise-corrupted measurements. 

In the model chosen for the spacecraft, only two basic defor- 
mation modes were considered in order to reduce the number of 
parameters to be identified. Test runs showed this hypothesis 
was quite acceptable in presence of damping. The filtering al- 
gorithm produced excellent results for a model which was not 
all that well suited, but requires great care in the choice of the 
initial values of the various state vectors and probability ma- 
trices. 


We may add that, apart from the reduction to two modes, the 
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(813) exact 


(85I,) exact 


(Sl sym ) exact 


Identification of 3 parameters for GEOS 


model is quite general : in a practical case, a preliminary sen- 
sitivity analysis may prove some of the parameters to be prac- 
tically constant. These parameters may then be discarded from 
the state-vector, thereby reducing the final order of the system. 

Moreover, recent results obtained for a satellite whose dy- 
namical model is very accurate give a further proof of the suit- 
ability of Kalman filtering for the type of probl:m considered. 

The major drawback of the algorithm is its time consumption, 
arising mainly from the integration of the error probability co- 
variance differential equation. A typical computer time for one 
filtering step is 10 ... 15 sec which practically excludes real-time 
operation. The method however proves useful for post-proc- 
essing of accumulated data. An alternate approach would be 
to implement a suboptimal filter, doing away with the integra- 
tions, which would probably be less accurate, but also less time 
consuming. 
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to further op:imize his payoff at the maximizing player's expense? 
would be the greatest lower bound for the maximizing player's payoff? To answer these 
questions, necessary conditions for parameter optimization for linear quadratic dif- 
ferential games will be derived. Iterative numerical techniques for determini.g optimal 


what 


Hence 


parameters will be developed. Search techniques which will locate a ‘‘small’’ region of 
uncertainty in which the optimal parameter must lie will also be discussed. 


Introduction 


A faudamental result of optimal control theory is that the 
optimal control of linear systems leads to an optimal performance 
cost that is a function of the system parameters. Therefore it is 
reasonable for the system designer to further minimize the per- 
formance measure by manipulating the values of certain param- 
eters if this can be done. This is what parameter optimization 
deals with. Above all, we could show that the above technique 
could be extended to linear quadratic zero sum differential 
games because the analytic methods derived for them are actually 
extensions of technique used in control, 

In the past years, several papers had dealt with this subject 
in system control. Brooks [1]? and Spang [2] used random search. 
Through that technique one can locate global optimal param- 
and therefore some- 


eters of nonlinear system. It is inefficient 


what impractical. IXokotovie and Heller [3] derived and com- 


pared two analytic expressions for V.J0, where w is the sys- 
Vial the cor- 
responding gradient of the criterion function, J. Subsequently, 
Unbehaven [4] applied gradient methods with one of Heller’s 


expressions to a system without control. 


tem parameter vector to be optimized and 


Hofer and Sagirow [5] 
and Boltyanskii [8] pointed out that Pontryagin’s [6] necessary 
Hofer 
also developed results for systems for which the parameters to 


condition for optimal parameter selection was erroneous. 


be optimized are not the coefficients of the state variables, 2x(/). 
Ahmed and Georganas [7] showed that Boltyanskii’s [8] results 
could easily be derived from Gamkrelidze’s general maximum 
principle {11}. After that, Georganas [9] presented imbedding 

IThis research was supported in part under AFOSR Research Grant 76- 
2958. 


2Numbers in brackets designate References at end of paper. 
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techniques for estimating optimal parameters, but he assumed 
the optimal parameter could be expressed analytically. To sum 
up, most of the above papers handled systems without control 
or systems with specified structure. This paper intends to gen- 
eralize the above techniques and extend them to linear quadratic 
differential games. 

We assume the minimizing player, U, has the opportunity to 
choose the values of certain parameters before the game starts. 
Under that what values of those 


circumstance, parameters 


should be chosen, if he wants to further optimize his own gain 
> 


at the maximizing player’s, V, expense? On the other hand, it is a 


“worst-case’”’ analysis problem for V. If he only knows some pa- 


rameters lie within a certain range, w; < w < w,, before the 
game starts, what would be the greatest lower bound for his ex- 
pected payoff? We first use the calculus of variations to derive 
the necessary conditions for the optimal parameter vector. Then, 
we develop the gradient, invariant imbedding, and quasilineariza- 
tion methods for determining optimal system parameters. After 
that, we discuss the searching technique for locating the initial 


guess region. 


Formulation of the Problem 


Consider a linear quadratic pursuit-evasion differential game 
with constant parameters: 
¢(t) = Ax(t) 4 


Gyu(t) Ga £): tlle) = Ze (1) 


where A is a constant parameter n X n matrix, x(¢) is the state 
n-vector, and v(t) and u(t) are “‘maximizing’’ control p-vector 
and ‘minimizing’ control m-vector respectively. The criterion 
function, J, is 


— vT(t)Ra t)|dt 


a = l 
a7™(T)M2x(T) + 
) 9 
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where M is an n X n semipositive definite symmetric matrix, 
and R; and R2 are m X mand p X p positive definite symmetric 
matrices, respectively. 

Necessary and sufficient conditions for the existence of a sad- 
dle-point are well-known [10]. The optimal payoff is 


1 
J(u*, v*) = - xoK (to) x0 (3) 


where the optimal strategy is: 
u*(t) = — Rr G,TK(t)z(t) 
v*(t) = R21G,TK(t)z(t) 


where 
p(t) = K(t)zx(t) 


and 

K(t) + K(t)A + ATK(t) + K(t)[(G:R2G,7 — G,R:G,7) 
K(t) = 0; 
K(T)=M_ (7) 


It is obvious that V has to use (5) regardless of the values of 
w (w, = ai;, the r elements of A which can be optimized) if he 
wants to achieve his goal. Similarly, U has to use (4). Conse- 
quently, U assumes that both (V and himself) will use the optimal 
strategy when he chooses the values for w before the game 
starts. Assume the existence of the saddle-point for ali admissible 
w, w will only change the value of the payoff. 


So 


J(u*, »*, w*) = min J(u*, v*, w) < J(u*, v*, w) (8) 
weRT 


Therefore, as long as both players know what w* is when the 
game starts, the optimal payoff to both players will be J(u*, 
v*, w*) if both players play optimally. 


On the other hand, even though V is not sure about how U 
picks the parameter w, if he uses the maximum policy, he will 
choose to maximize his cost against the worst possible strategy 
which U That is, assuming U’ pick u*(t) and w%, 
he will use r*(/) and w*. Then his expected payoff is also J(u*, 


»* 


will choose. 


, w*), the greatest lower bound. 
It is obvious that we can use the augmented state approach to 
find the optimal constant parameter vector w* by considering 
the r elements of w as additional state variables which are free 
at both ends. 
Variational Approach [14] 

Consider the system defined by (1) and (2). The Hamiltonian 
is 


l 
H(2x(t), u(t), v(t), p(t), w, t) = . [u?(t)Riu(t) — v7(t)Rv(t)] 


+ pT(t)[Ax(t) + Giu(t) + Gw(t)] + AMO) = (9) 
By the calculus of variations 
T 


éJ = [Mz(T) — p(T)|6rr + AT(t)Sw 


to 


+ p(t) + or + | z(t) — 
J to Ox 
‘ 0H jt 
+ | io + : | bw 
Ow 
OH 
= Ou 


oH }F 


ap | °P 
b(t) ay bX 

~(t) — ) 

+ | we) ad 

T T 

| bu + | | ao 


(10) 
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If r(t)eR*, u(t)eR™, v(t)eR?, weR’, then the necessary condition 


for optimality is 6J = 0. So after some manipulation, we get 


‘ OH 
A(t) oo aw X(to) = (0 and X(T) = 0 (11) 


w(t) = 0 
Z(t) = Az(t) + (G2:R27G,7 — GiR,"G,")p(t); 


(12) 

r(to) = % (13) 
and 

p(t) = — ATp(t); 


It can easily be shown that 


p(T) = M2(T) (14) 


h(t) = — pilt)ax;(t); Ane) = 0 and A(T) =0 


where w; = a1; for k = 1, 


(15) 


Tr. 


X(t), the Lagrange multiplier, is used to adjoin (12), the dif- 
ferential equation describing the unknown parameter, to the 
criterion function. 

Assume the existence of the optimal strategies u*(t) and 
v*(t) for all admissible w, if we let dw be constant throughout the 
interval of integration, for all te[éo, 7'], we get 

T OH 


T 
6J = [dw] i dt 


= 2 bux f ms 
[pr(t)x;(t)|xdt = 0 


(16) 
k=l 

Therefore (15) or (16) is the additional necessary condition 
for the optimal parameter vector. However, if w is allowed to 
vary not in the entire space R’, but only in some closed domain 
[w:, w.], having a piecewise smooth boundary, then the addi- 
tional necessary condition for w to be the admissible optimal 


parameter is 
rToH |t 
se ee eae 
t | Ow 


where 6% must be in the feasible direction. Subsequently, we 
will develop several numerical techniques to find the global op- 
timal parameter vector if w is unbounded or the admissible 
optimal parameter vector if w is bounded. 


(17) 


The Modified Gradient Method 
By inspecting (16), we could define 
= — Tdi! 


dws = wit! — wt 


where 


T 
bJ,f = f. [pi(t)a;(t)Ja*dt 
to 


SJ. = > 6w,bJ 5 = — ) 3 T(J 5)? 


k~l k=l 


(20) 


(20) indicates that the gradient method is always searching in a 
feasible direction. We stop updating wy if léJ;‘ll < € or if we 
reaches the boundary after some iterations. 
Invariant Imbedding [15] 
Let 
[z7(T)w7(T))? = M(S, T) = N(T) — G(T)S (21) 


where 
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S = [p"(T)A"(T)]* 
We use the invariant imbedding equation 


aM 
aT 


OM 


—— 9(] =f 
+ 36 g(M, 8, T) = f(M, 8, T) 


[é7(t)wT(t)|7 = f{(x7, w?)?, (p™, AT)?, (24) 
[pT (t)AT(L)IT 


Using (12) to (15), we substitute (24) and (25) into (23) and get 


g{(z™, wT), (pt, AT)*, a (25) 


Nit) — GS — GWPWS = FING — EWS) 


+ [B,R21B,7 — ByR,1B,7|S (26) 


4 ‘0 
--i--|(n+r)X (n+r) 
0 :@ 


G2 
B, = --[e +r)xXm B, = | .... In +r) X p 


0 0 


The m; corresponds to the z; defined in (15). Ignoring terms 
higher than first order in S, we get an n + r vector differential 
equation, N(¢), for terms without S and an n + r square matrix 
differential equation, G(t), for the terms with S. The initial con- 
ditions are N(éo) = [ao7C7]’ where C is the initial or updated 
value for the unknown parameter w and 


0 id 
G(to) = | -- 1 -- 
dT \d 


where G(to) is a symmetric matrix with an n X n zero matrix and 
weighting matrices d and d. Though the elements of d and d are 
arbitrary constants, values which are too high may result in di- 
vergence, whereas, values which are too small may cause very 
slow convergence. 
By partitioning 
nXn | 
Gi(T) | G(T) 
ene 7) ates ccs ft cane 
G(T) | G(T) 
rXn 


Ni(T) nX1 


and using p(7') = Mz(7') and \(7’) = 0, we can verify that 


a(T) = U + G(T)M}N\(T) (27) 


w = NT) — G2(T)M2(T) (28) 

w obtained in (28) will be used as the initial condition for C in the 

next iteration. The stopping criterion is max ||\;‘(7')ll < « We 
k 

can reduce computations if we use max |lw,'+! — w;'ll < 6 as the 

k 
stopping criterion. By so doing, we do not have to find the tra- 
jectories x(t), p(t) and A(t) at each iteration. 


Quasilinearization [12] 


Employing a first order Taylor series expansion for (12), (13), 
(14), and (15), we get the (¢ + 1)th estimate for 8 from the ith 
estimate 


Biss = TiBign + Qi (29) 
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v(t) AAr()) | 


Ow 
mt O[Az(t)] ; 
Ow . 
[pr (taste 
L(t) ji L 0 J 
_ AA) 7 
Ow 
(G2R2G,7 - G,R,14,") A 0 ote 
—2xj(t)e —pi(t), 0 0 
0 0 O 0 


x(t) 
A(t) 


where Bix; = 














—At 0 O 








— 


Therefore, the homogeneous equation is 
Bin® = TiBin® 


and the n + r sets of initial conditions are 


[B¥1(to)B¥2(to), ..-, B¥nsr(to)] = 0 ~2(n+r) X (n+ r) 


where the first n X n and the last r X r identity matrix are values 
for p#(to) and w* (to), respectively. 
The nonhomogeneous equation for the particular solution is 


Bis? = TBin? + 4 (31) 
where the only set of initial conditions is 
B?(to.) = [072970707]? 
And it can be shown that 


p¥(T) — Mz*(T) | | p¥ase(T) — Mr®ay(T) | 
Ci = pee 
A#(T) held Mase(T) i 


Mz?(T) — p?(T)~ 
coms AP( T) , 


pi*l(%) = [Ci, .. 


2(n + r) vector (32) 


(33) 


where 
"> Cala r 


and 


wit = (Casa, ces Casrliga™ 


(34), (35), and 2 give us sufficient information to solve the orig- 
inal problem directly. If || (¢)ll < €, we get the optimal solution, 
otherwise, we use them as the nominal trajectory for the next 
iteration. Since, [’ is a sparse matrix we can often reduce the 
dimension by partitioning it into several subsystem equations. 


Searching for the Convergent Region 


Since the augmented state approach, invariant imbedding and 
quasilinearization are all locally convergent techniques, it is 
essential to determine a local region of uncertainty in which the 
optimal parameter vector must lie. Here, we only discuss a 
scalar case to illustrate our techniques. It is well-known that a 
linear quadratic performance measure is convex. Therefore, the 
three-point pattern which is shown in Fig. 1 can give us an idea 
of where the region of uncertainty is. The optimal parameter is 
to the left, to the right, or in between those points. And it is 
worthwhile to point out that we only need (3) and (7) to find 
J(w). Once the region is located in between two points, we can 
further reduce the region by applying Fibonacci search [13], and 
it is obvious that we could generalize the above technique to 
vector case. 
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Fig.1 Three-point pattern 


Example 


Given: System Equation 
‘pone ss ie w ] a(t) | Pe Jue 
t2(t) 0.5 0.5j] x(t) 1.5 
‘ 1 he a(t) | _ [2.0 
0 X2(to) 2.0 


and the performance measure 


(36) 


1 
J = - a7(T)M2(T) + 


- 


1 fos 
f, [u2(t) — v2(t)|dt (37) 
2 Jo 


where M is an identity matrix. 


The gradient method 
9 
—f f. [pi(t)ae(t)|dt 
0 


Using wo = — 0.5 where J = — 1.338 
where J = 1.484 in 10 iterations. Using w = — 0.1, we get 
w* = — 1,338 in 11 iterations. The 7 is determined by 7 = 0.3 
+ 0.05 N where N is the iterative number. 

Invariant Imbedding 
Since n3(¢) = 0 and n;(t) = w, we get 


ined _ ‘te nN3 m(t) i; ciel re ° 
fio(t) 0.5 0.5 ]] no(t) | | nolte) 2.0 


and 

gu(t)gie(t)gis(t) —1 
ga (t)q20(t)ges(t) w 
gsi (t)gso(t)gsa(t) no(t) 


bw = (38) 


1.867, we get w* = 


Jur(t)gro(t)gis(t) 
Gar (t)Goo(t )gos(t) = 
Gar (t)@so(t)gs3(t) 


gai(t)ne(t) 9s2(t)ne(t) 9s3(t)no(t) —] 
+ 0 0 0 + | 0.5 
0 0 0 


gu(t) 
galt) 
gsi(t) 


gio(t) 
922(t) 
gso(t) 


gis(t) 
J23(t) 
9s3(l) 


d 


G(to) = 0 0 
dda 


(40) 


The integer d for 1 < d < 8 has been used as the initial condition. 
Using w) = — 0.5 where J = 1.867, we get w* = — 1.345 in 11 
iterations for d = 4.0. We need 15 iterations for d = 1.0, and 
the value for w fluctuates around the optimal value when d > 7.0. 
The stopping criterion used is |lw**+! — will < 0.002. 
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Quasilinearization 


Since w*(t) = w‘(te) and all equations are independent of X(t), 
therefore, with numerical values, (29) can be rearranged into 


pr(t) 1.0 0 
P(t) ) “t 0 
Z(t) ; w 
A(t) Jia —2,(t) —pr(t) ji Laa(t) 


0 0 
— p(t) pi(t)w 
+ a(t) witl + | —wa(t) 
0 0 
0 F pr(t)xe(t) Js 


pit) 
pr(t) 


Using wo = — 0.5 as an initial condition, we get the optimal 
value w* = — 1.338 in 6 iterations. 

All computations were performed on UCLA’s IBM 360/91 
computer, and all integrations used a fourth-order Runge-Kutta 
method with At= 0.005, 


Discussion 


Several models have been used in each technique, and all 
optimal parameters settle on a “stable value” unless the system 
is unstable for all admissible parameters. Also, it can be easily 
shown that the necessary condition (11) and the numerical 
techniques developed for differential games can be applied to 
systems control problem with minor modification. In addition 
to always searching in the feasible direction, the gradient method 
initially converges rapidly ; however, it slows down as it gets close 
to the optimal value. So it may be used as a starting procedure 
and quasilinearization may then be utilized to close in on the 
solution. Above all, using the initial searching technique to 
locate the convergent region will ensure that both the invariant 
imbedding and quasilinearization method will converge to the 
true optimal parameter vector w*. Furthermore, it is worth- 
while to point out that all three techniques have similar stopping 
criterion since ||\,‘(7')I| = WéJ;‘ll. Intuitively, the developed 
techniques can be extended to time varying parameter, w(t), 
where differential equation describing w(t) will be used instead 
of (12), and (11) will be the additional necessary condition. 

Anyway, this paper only solves half of the parameter optimiza- 
tion problem in differential games. How the maximizing player 
chooses the values of certain parameters in order to achieve his 
goal. In this case, we should restrict the range of the values of 
the parameters to ensure the game problem is meaningful, that is, 
the game can be terminated. 
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Dynamics of an Electrohydraulic Stepping Motor 





M. J. Vilenius' 


In this paper a nonlinear mathematical model is presented for 
a modified commercially constructed electrohydraulic stepping 
motor. With the aid of the model, the dynamic performance of 
the motor is examined and the quality of the model tested. 


Nomenclature 


A; = area of servovalve spool 
Ay = area of nozzle 
saturated value of error amplifier voltage 
bulk modulus 
viscous damping coefficient of servovalve armature 
viscous damping coefficient of servovalve spool 
viscous damping coefficient of hydraulic motor and 
load 
viscous damping coefficient of electric stepping motor 
leakage coefficient 
volumetric displacement of hydraulic motor 
diameter of servovalve spool 
frequency of clock pulses 
total inertia of servovalve armature 
total inertia of hydraulic motor and load 
total inertia of electric stepping motor 
spring constant 
torsion spring constant 
flow-pressure coeilicient of nozzle valve 
gain of error amplifier 
coefficient of potentiometer 
flow gain of nozzle valve 
gain of servo amplifier 
torque constant of torque motor 
coefficient associated with Bernoulli force 
total mass of servovalve spool 
maximum torque of electric stepping motor 
number of rotor poles 
supply pressure 
length 
time 


time between two clock pulses 


input signal 

coefficient associated with orifice area 

half of total volume of fluid under compression in 
hydraulic motor 

contained volume at each end of servovalve spool 


1Act. Associate Professor, Dept of Mech. Eng., Tampere University of 
Technology, Tampere, Finland. 
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angular deflection of hydraulic motor shaft 

load pressure difference 

displacement of servovalve spool 

pressure difference between ends of servovalve spool 
angular deflection of servovalve armature 

n Qs 

discharge coefficient 

density of oil 

angular deflection of electric stepping motor shaft 


Introduction 


This paper is based on research, which was done in connection 
with designing the control system for the table of a special drill- 
ing machine. Because the control signal was digital and the 
inertia load rather large, an electrohydraulic stepping motor was 
one of the most viable and economical alternatives. In addition, 
the control system had to be quite adjustable which meant that 
a stepping motor construction with an electrohydraulic position 
control servo seemed to be a better alternative than a construc- 
tion with a hydromechanical position control servo. -However, 
for the present purpose a commercial stepping motor with an 
electrohydraulic position control servo was too inaccurate, and 
its steady state angular velocity was unnecessarily high. There- 
fore the commercial unit was modified to fulfill the requirements 
of the present application. It was also important to be able to 
predict the effects of different factors on the dynamic perform- 
ance of the stepping motor and therefore a rather accurate math- 
ematical model was developed to assist in designing the special 
drilling machine. In the model the most important nonlinearities 
have been taken into account. The dynamics of the stepping 
motor are examined with the aid of this model. The quality of 
the mathematical model is tested by comparing it with the ex- 
perimental results. 


The System 


A schematic diagram of the stepping motor construction is 
shown in Fig. 1. The digital signal given by the control unit is 
changed to an analog signal (potentiometer shaft angle) with the 
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Fig.1 Schematic diagram of system 


MARCH 1977 / 63 





aid of control logic and the electric stepping motor. This signal 
is amplified by the electrohydraulic position control servo. The 
output signal of the entire system is the angular deflection of the 
hydraulic motor shaft. 

The commercial electrohydraulic stepping motor had insuf- 
ficient accuracy for the abovementioned application. The 
greatest error was due to the threshold region of the two-stage 
electrohydraulic servovalve. Since the steady state angular 
velocity was unnecessarily high, it was possible to replace the 
servovalve with a smaller unit which still fulfilled the steady 
state angular velocity requirements. The new servovalve had a 
smaller threshold region and also had better dynamic charac- 
teristics, so that it was possible to increase the loop gain of the 
position control servo in the modified stepping motor. 

Because of the modification, the accuracy of the unloaded 
electrohydraulic stepping motor became about five times better 
than that of the original commercial motor. In addition, re- 
sponse times were found to decrease slightly, a consequence of 
the improved dynamic characteristics of the new servovalve. 


The Mathematical Model 


Nonlinear differential equations describing the dynamics of 
the electrohydraulic stepping motor have previously been de- 
veloped [1, 2, 3, 4, 5).2. The nonlinearities taken into account 
were the torque of the electric stepping motor, the voltage limit- 
ing of the error amplifier, the Bernoulli force acting on the 
servovalve spool and the load flow through the servovalve. From 
the following set of nonlinear differential equations, the system 
equations in the state-variable form were solved. 

The state function was found to be 


x(t) = F(x(t), w(t)) 
where 
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when |(1/n)z» — 2| < a (the operating region of the potentiom- 
eter) and —a < K.K,((1/n)a_ — 21) < +a (the linear operating 
region of the error amplifier) 


*Numbers in brackets designate References at end of paper. 
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(a + 1)(w/4), whenaT <t < (a+ 1)T anda = 0,1, 
2, .... In the area |(1/n)ay — m| < mw and |K-K,((1/n)2» 
— 2)| > a (the voltage limiting region of the error amplifier) 
the 8th row of the state function (1) changes into form 


a-s)] 


but in the other parts the state function remains unchanged. 
The output expression is a scalar equation 


y= %. (3) 


Dynamic Performance 


Based primarily on the technical characteristics given by the 
manufacturers, the fixed parameters of the mathematical model 
(Table 1) have been ascertained. Using these fixed parameters, 
the way in which some of the adjustable parameters affect the 
dynamics of the electrohydraulic stepping motor has been ex- 
amined. In the mathematical model the output signal (i.e. the 
angular deflection of the hydraulic motor shaft) is shown as the 
state variable 2). All responses have been calculated by using 
the simulation language MIMIC of the computer UNIVAC 
1108. 

The effect of the gain K, on the response of the state variable 
Z; was examined with the hydraulic motor unloaded (J, = 
0.00293 kgm? = the inertia of the rotor of the hydraulic motor), 
having the input signal of one clock pulse. The responses of the 
state variable x, to one clock pulse are shown in Fig. 2 for three 


Table 1 Fixed parameters of the mathematical model 
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= 500 Ns/m 
(motor without load) 
(motor with load) 
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Fig. 2. Single pulse responses of x:, showing the effect of gain K, 
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different values of K,. As expected the response time of the state 
variable x, decreases while the gain K, increases. The oscillation 
of the responses is caused by the oscillation of the electric stepping 
motor shaft. When the gain K, decreases the total gain of the 
position control servo decreases, and the effect of the oscillation 
of the electric stepping motor shaft on the response of the state 
variable 2; is damped. 

The effect of the frequency of the clock pulses on the response 
of the state variable x: was examined with the gain K, = 1.5 


Journal of Dynamic Systems, Measurement, and Control 


when the hydraulic motor was unloaded. The responses to three 
clock pulses at the frequencies 15 Hz and 200 Hz are shown in 
Fig. 3. It can be noticed that at certain frequencies (e.q. f = 200 
Hz) the response doesn’t oscillate but rather is optimally damped. 

The effect of the inertia load on the response of the state 
variable 2; was examined with an input signal of one clock pulse 
and the gain K, = 0.3. For the responses to one pulse (Fig. 4) 
the hydraulic motor was either unloaded (J,, = 0.00293 kgm?) 
or the inertia load was relatively high (J,, = 1.23 kgm?). A high 
inertia load caused an increase in the response time and an over- 
shooting of the response. 

In order to prove the quality of the mathematical model the 
experimental responses of the state variable x, are also shown in 
the figures. Slight differences appear between the theoretical 
and experimental responses, but when designing the control 
systems in practice the mathematical model can be regarded as 
completely sufficient in all cases. 


Conclusions 


A nonlinear mathematical model is presented for an electro- 
hydraulic stepping motor. With the aid of the model, the way 
in which some of the adjustable parameters affect the dynamics 
of the motor has been examined. 

The frequency of the clock pulses affect strongly the damping 
of the electrohydraulic stepping motor. At certain frequencies 
the response doesn’t oscillate but rather is optimally dampened. 
The increase of the total gain of the position control servo in- 
creases the oscillation of the response and decreases the response 
time. The high inertia load causes the increase in the response 
time and the overshooting of the response. 
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“Cybernetic Analysis of Societal Systems and World Modelling, Part III,” 
P. N. Rastogi, Ibid., Vol. 27, No. 4, Oct. 1976, pp. 136-141. 

“Modeling and Optimization of Economic Systems in Socialist Countries," 
K. Cichocki, J. of Dynamic Systems, Measurement, and Control, (USA), Trans. 
ASME, Vol. 98, Series G, No. 3, Sept. 1976, pp. 252-259. 


Part II," 


MANAGEMENT 


“Periodic Job Scheduling in a Distributed Processor System,"’ M. J. Gonzalez 
and J. W. Soh (USA), IEEE Trans. on Aerospace and Electronic Systems,"’ 
Vol. AES-12, No. 5, Sept. 1976, pp. 530-535. 

“Policy Analysis Modeling for Community Development in Small Cities,"’ 
F. Dicesare, D. A. Ensminger and K. R. Carlisle (USA), JEEE Trans. on 
Systems, Man and Cybernetics, Vol. SMC-6, No. 8, Aug. 1976, pp. 515-531. 

“An Experimental Method of Determination of Optimal Maintenance 
Schedules in Power Systems Using the Branck-and-Bound Technique,” G. T. 
Egan, T. 8. Dillon; and K. Morsztyn, Ibid., pp. 538-547. 

“Scheduling Parallel Production Lines with Changeover Costs: Practical 
Application of a Quadratic Assignment/LP Approach,’’ A. M. Geoffrion and 
G. W. Graves, Operations Research, Vol. 24, No. 4, July /Aug. 1976, pp. 595-610. 

“Analysis of Single Critical Number Ordering Policies for Perishable In- 
ventories,"" M. A. Cohen, Ibid., pp. 726-741. 

“Energy, Employment and Growth,"’ (in Swedish), A. Bigsten and L. 
Hjelmarsson, Ekonomisk debatt (Sweden), Vol. 4, No. 3, 1976, pp. 195-206. 

“‘Model for Maintenance Work in a Socialist Economy" (Yugoslavia; in 
Swedish), E. Rejec and F. Mencinger, I ndustriell Teknik (Sweden), Oct. 1976, 
pp. 30-33. 


ACOUSTICS AND NOISE 


“Characterisation by Energy of Short Duration Noise in Logical Systems,” 
J. L. Aucouturier, J. P. Dom, G. Lacroix, J. J. Riviere, L’’Onde Electrique 
(France), Vol. 56, No. 10, Oct. 1976, pp. 397-402. 

“Noise Measurements for Very Low Impedances,"’ G. Lecoy, D. Rigaud, F. 
Segui, Ibid., pp. 403-406. 

“Sound Analysis by Mini-Computer in Conversational Mode,” 8. Chabrel, 
3. Charbonneau, L’Onde Electrique (France), Vol. 56, No. 8/9, Aug./Sept., 
976, pp. 358-366. 

“A Model of Community Noise Pollution,’’ R. E. Burke, (USA), Simulation, 
Vol. 27, No. 3, Sept. 1976, pp. 87-96. 

“Speaker-Recognition Using Orthognal Linear Prediction."’ M. R. Sambur 
(USA), IEEE Trans. on Acoustics, Speech and Signal Processing, Vol. ASSP-24, 
No. 4, Aug. 1976, pp. 283-289. 

“Signal Analysis by Homomorphic Prediction,’’ A. V. Oppenheim and J. M. 
Tribolet, Ibid., pp. 327-332. 

“‘A Comparative Performance Study of Several Pitch Detection Algorithms,” 
L. R. Robiner and C. A. McKumegal, Ibid., Vol. ASSP-24, No. 5, pp. 399-418. 

“Maximum Likelihood Pitch Estimation,” J. D. Wise and T. W. Parks, 
Ibid., pp. 418-423. 
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book reviews 


APPLIED MECHANICS: MORE DYNAMICS, 
John Wiley & Sons, Inc., 1976, 244 pp. 


by Charles E. Smith, 


REVIEWED BY GEORGE R. SPALDING! 


When this book arrived from the publisher, the reviewer was 
planning a one quarter course in advanced dynamics which was to 
contain applications to inertial systems. The intention was to 
present rigid body dynamics, rotating coordinate systems, linear 
coordinate transformations, Lagrange’s equations and gyroscopic 
motion, and then to apply these concepts to gyroscopic devices 
and inertial systems. The applications were to be taught from 
selected technical articles. The major problem, until Professor 
Smith’s book arrived, was finding a text that covered these basic 
topics in a continuous and concise manner. 

It contains 
five chapters, the first four of which cover the material necessary 


More Dynamics met the requirements perfectly. 


for an understanding of gyroscopic systems. 

The first chapter starts with the time derivative of a vector in 
a rotating space, develops the velocity and acceleration rela- 
tionships, briefly discusses particle kinetics, and ends with gen- 
eral motion of a rigid body. Rotating frames of reference are 
always emphasized. 

The next chapter begins with basic linear algebra and goes 
through coordinate transformations and the determination of 
principal directions. For most engineering students, this material 
will not be new; however the examples and the problems at the 
end of the chapter serve as an excellent review. 

The third chapter contains a fine discussion of rigid body 
rotational displacements, Eulerian angles and small angle ap- 
proximations. This chapter also presents angular momentum, 
kinetic energy, Newton’s second law, and ends with a limited 
discussion of the top problem. 

The development of Lagrange’s equations, from virtual work, 
is treated in the fourth chapter. Degrees-of-freedom, generalized 
coordinates, potential functions and the derivation of the equa- 
tions are clearly and concisely presented. This is followed by a 
brief discussion of the Hamiltonian, energy integrals, and the 
top problem, this time approached from the Lagrangian. 

The book’s final chapter is an excellent undergraduate-level 
discussion of the dynamics of engineering systems. The important 
The 


emphasis is on linear constant parameter systems and _ the 


points are here, presented in a straightforward manner. 


meaning and consequences of linearity. The transition from the 
discrete model to the spatially distributed one is very well made. 
There is a discussion of linearization, and the chapter concludes 
with a few pages on some common nonlinear characteristics. 

The first four chapters provide an excellent background for 
further work on inertial navigation. The notation used to dis- 
tinguish reference frames was particularly well chosen. For ex- 
ample, if a and 6 designate two coordinate systems rotating with 
respect to each other, then the equation 


a B 
A=A+,_2sXA 


‘Department of Engineering, Wright State University, Dayton, Ohio. 
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indicates that the time derivative of A in a-space equals the time 
derivative in 8-space plus the angular velocity of 8 as observed 
from a, crossed with A. The students liked this notation very 
much, 

While the final chapter is excellently done, it does not fit 
naturally at the end of the other four. It does not build on them, 
nor apply the concepts developed in them. A chapter which did 
would, in the reviewer’s opinion, improve this otherwise excellent 
text. 


NOTES ON NONLINEAR SYSTEMS, by J. K. Aggarwal, Van No- 
strand Reinhold Notes on System Sciences, 1972, 201 pp., 
Paperback. 


REVIEWED BY S. H. JOHNSON! 


As part of a series of such books, Van Nostrand Reinhold has 
published the classnotes for a one-semester interdisciplinary 
course on the methods of qualitiative quantitative and numerical 
analysis of systems of nonlinear ordinary differential equations. 
These notes have been used with first-year graduate students 
and seniors at the University of Texas. It has been used for a 
half semester course segment for Lehigh graduate students. The 
book is intended to complement the presentation in two other 
books in the same Van Nostrand-Reinhold series: Notes for a 
First Course on Linear Systems, by Polak and Wong and Notes 


for a Second Course on Linear Systems by C. A. Desoer. 


The book contains the things one would expect: Phase-Plane 
Methods, Stability of Linear and Nonlinear Systems, Limit 
Cycles and Ultimate Bounds, Approximation Methods and 
Numerical Integration Algorithms. The presentation of the ma- 
terial is particularly appealing. Many topics are developed in 
same fashion as they evolved historically. Examples and -exer- 
cises are taken from the literature with the sources cited. Un- 
fortunately this technique has not been well executed in detail. 
Information is omitted from the exercises which would greatly 
reduce the drudgery. The reader is asked to sketch trajectories 
for Rapoports’ arms-race equations without being provided 
numerical values for the parameters. If the reader consults the 
cited reference he fails to find numerical values which would 
make an interesting and demonstrative exercise and fails to find 
the equations attributed, albeit somewhat obliquely, by the 
author to that reference. 

There are a few mistakes or omissions in the book but too few 
to be annoying. All of the numerical results in the computer 
methods chapter have been reproduced without discrepancy. The 
book is physically ill-suited to use as a text. Pages fall out so 
often that students resort to drilling holes in the book and con- 
verting it to loose-leaf form. The Editor’s Note inside the back 
covers says that a clothbound library edition is also available. 

Despite the shortcomings, the book is useful and will be used 
again by the reviewer. 


1Associate Professor of M.E. and Mechanics, Lehigh University, Bethlehem, 
Pa. 
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1977 JOINT AUTOMATIC CONTROL CONFERENCE 
Hyatt Regency Hotel 


San Francisco 
June 22-24, 1977 


The theme of the 1977 JACC is ‘Control Theory and Applications in the Ser- 
vice of Local Industry,” and includes tutorial state-of-the-art and general sessions 
in the areas of System Modeling and Identification, Optimal Self-Organizing, 
Learning and Hierarchical, Decision and Control, Pattern Recognition and 
Linguisitic Methods applied to Large Scale Systems such as Economic, Trans- 
portation, Power, Seismic Data and Oil - exploration as well as Space and Air- 
craft Reliability, Health Care and Bioengineering Systems, Robotics and Indus- 
trial Manipulators and the use of Microprocessors in Control and other related 
areas. 


Call for Papers—1977 WAM 


The Automatic Control Division of ASME is planning sessions for the 1977 Winter Annual 
Meeting, to be held in Atlanta, Georgia, November 27-December 2, 1977. Selection of topics 
for several sessions is stimulated by the results of recent Survey of ACD membership which 
indicated strong cross links between control and areas represented by the Design, Management, 
and Process Engineering Divisions. Sessions planned on transportation and energy reflect the 
high level of current research activity that benefit from advances in control theory and applica- 
tion. 

Papers are requested in the following areas: design automation and optimization, robotics and 
manipulators, computer control of industrial processes, productivity, stability and control in rail 
transportation, and control and optimization of energy supply, conversion, and end-use technol- 
ogy. Case studies of application of results in the above areas in industrial situations are welcome. 
Sessions will also be arranged for timely contributed papers not included in the above sessions. 
Send five copies of papers to Dr. Karl N. Reid, Sr. Technical Editor, JDSMC School of Mechanical 
& Aerospace Engineering, Okla. State University, Stillwater, OK 74074, 405-624-5900 by April 
1, 1977. 


Fluid Control Systems Committee—Call for Papers 


The Fluidics Committee has changed its title to “Fluid Control Systems Committee” and has 
expanded its interests to include theoretical and system design aspects of all fluid control systems. 
A principal interest of the committee will still be in fluidies; however, moving-part fluidic control 
systems using hydraulics or pneumatics in either digital or analog applications are now included. 

In keeping with its expanded role, the Fluid Control Systems Committee is soliciting papers in 
the area of: component development or description, systems design, control applications, or 
modeling and simulation of fluidic, pneumatic, or hydraulic systems. The theme for ’77 WAM 
is noise, thus papers concerning signal noise characteristics or analysis, or airborne acoustic noise 
are of special interest. 

Authors should mail five copies of their manuscript before April 1, 1977 to Dr. Robert L. Woods, 
Associate Editor JOSMC, JFE, Mechanical Engineering Department, University of Texas at 
Arlington, Arlington, TX 76019. 
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Please 

Submit manuscripts 
for review and 
address all 
communications 

to: 


Prof. Karl N. Reid 

School of Mech. & Aero. Eng. 
Oklahoma State University 
Stillwater, Okla. 74074 


information for authors 


«The purpose of the Society is to disseminate technical information of perma- 
nent interest resulting in the publication of a series of Transactions quarterlies. 

«Papers of permanent interest having been selected for publication in the ASME 
Transactions are subject to a Page Charge of $35 per page, which will be invoiced, 
through the author, to the author’s company, institution, or agency at the time 
page proofs are submitted to the author. Publication is not dependent upon pay- 
ment of the page charge. 

-Manuscripts should be submitted in final form to the Editor. Each manuscript 
must be accompanied by a statement that it has not been published elsewhere 
nor has it been submitted for publication elsewhere. A paper which would occupy 
more than 6 pages of the Journal will be returned to the author for abridgment. 

«The author should state his business connection, the title of his position, and 
his mailing address. A short abstract (50 to 100 words) should be included on 
the first page immediately preceding the introductory paragraph of the paper. 

Five copies of the manuscript are required. One of these must be a carefully 
prepared printer’s copy, typed in double spacing, on one side of the page only, 
with wide margins, on 22 by 28 cm (8)4 by ll-in.) opaque white paper. Copied 
manuscripts, if prepared with exceptional care, can be accepted provided they are 
completely edited. 

- All manuscripts submitted to ASME must use SI (Metric) Units in text, figures, 
or tables. In addition to SI units, English units may be included parenthetically. 

- All mathematical expressions should be typewritten. Greek letters and other 
symbols not available on the typewriter should be carefully inserted in ink. Care 
should be taken to distinguish between capital and lower-case letter, between 
zero (0) and the letter (O), between the numeral (1) and the letter (/), etc. A letter 
representing a vector cannot be printed with an arrow above or below it. The 
letter should be underscored with a single wavy line wherever it appears in the 
text, to designate boldface type. 

«A list of symbols carefully marked for the use of the editor (thus: x, Greek 
l.c. kappa), if it has not been included in the body of the paper, should accompany 
the manuscript on a separate page. 

Numbers that identify mathematival expressions should be enclosed in paren- 
theses. Numbers that identify references at the end of the paper should be en- 
closed in brackets. Care should be taken to arrange all tables and mathematical 
expressions in such a way that they will fit into a single column when set in type. 
Equations that might extend beyond the width of one column (fractions that 
should not be broken or long expressions enclosed in parentheses) should be re- 
phrased to go on two or more lines within column width. Fractional powers are 
preferred to root signs and should always be used in more elaborate formulas. 
The solidus should be used instead of the horizontal line for fractions whenever 
possible. 

«The text of an ASME paper normally should not exceed 6000 words (6 printed 
pages in a journal) or equivalent. A Technical Brief or Brief Note should not ex- 
ceed 1500 words or equivalent. Discussions should not exceed 500 words, and clo- 
sures normally should not exceed 250 words per discussion. In computing equiva- 
lence, a typical one-column figure or table is equal to 250 words. A one-line equa- 
tion is equal to 30 words. The use of a built-up fraction or an integral sign or sum- 
mation sign in a sentence will require additional space equal .to 10 words. 

-Originals and four copies of figures must accompany the manuscript. Line 
drawings should not be larger than 22 by 28 cm (814 X 11-in.)and should be planned 
for reduction to column width. Lettering should be large enough to be clearly 
legible when the illustration is reduced. The originals of line drawings must be in 
India ink or white or pale blue tracing paper or tracing cloth. Photographs of 
equipment or test specimens must be glossy prints and should be used sparingly. 
Captions for figures should be typed double-spaced and included as the last page 
of the manuscript. The figure number and author’s name should be written in 
the margin or on the back of each illustration. 

- Titles of papers should be brief. 

- Authors can obtain copies of the ASME Manual MS-4, ‘“*An ASME Paper,” from 
the Editor and are urged to do so before drafting their papers in final form. 

-An author is entitled to 25 preprints free of charge (in the case of two authors, 
15 each; three or more authors, 10 each). Larger quantities of preprints or re- 
prints can be ordered from the ASME Reprint Department, 345 East 47th Street, 
New York, N. Y. 10017. Quotations will be sent on request. 
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NEW BOOKS 


Just Published. . . 


Cavitation and Polyphase Flow Forum—1976 
Editor: Robert L. Waid 


The 17 papers which comprise this Forum present a wide range of primarily experimental studies of cavitation 
and polyphase flows. They include the effect of polymers on cavitation, boundary layer separation, and jet flows; 
modeling of cavitation on propellers and struts; measurement of interfacial waves; and impact and erosion of 
solids by liquids. The papers are arranged by subject and were presented at the 1976 Joint Gas Turbine and 
Fluids Engineering Conference held in New Orleans, Louisiana. 

1976. Book No. 100097, 56 pgs. 

$10.00, ASME members $5.00 


Centrifugal Compressor and Pump Stability, Stall and Surge 
Editors: P. C. Tramm and R. C. Dean, Jr. 


Gives practical considerations relative to centrifugal compressors and pumps, at least regarding their application to 
gas turbine engines. However, as the eleven papers of this Symposium demonstrate, there is today no fundamental 
understanding, nor an adequate theory, for predicting centrifugal turbomachine stability and flow breakdown. 
These papers were presented at the 1976 Joint Gas Turbine and Fluids Engineering Conference held in New 
Orleans, Louisiana. 

1976. Book No. 100098, 202 pgs. 

$20.00, ASME members $10.00 


Numerical/Laboratory Computer Methods in Fluid Mechanics 
Editor: A. A. Pouring 


The importance and significance of fluid mechanics has been known for centuries. Over the years, efforts focused 
mainly on analytical solutions of the governing differential equations, but as these became exhausted, more and 
more interest developed in “approximate” numerical methods of solution. 


The 24 papers which are included in this volume were presented at the 1976 ASME Winter Annual! Meeting. They 
cover navier stokes equations; boundary layer equations; heat transfer and boundary valve problems; machinery 
and transient flow; and laboratory methods. These papers review the current status of broadly selected areas of 
fluid mechanics and will benefit engineers working in this area by providing a permanent record. 

1976. Book No. G00106, 389 pgs. 

$35.00, ASME members $17.50 
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1976 Proceedings of the Fourth Intersociety Confer- 
ence on Transportation 


The papers contained in this volume place greater emphasis on the impact of regula- 
tion and legislation, energy and fuel, policy and planning, command and control, 
design and operation, and the total environment. They are important to a broad range 
of individuals and organizations — from planners and operators to users of transpor- 
tation systems. Some of the subject areas covered are: air cargo; highway and rail 
freight; passenger modes; water transportation; mass transit; pipelines; multi-modal 
energy transportation; economics and costs; policy and planning; command and con- 
trol; environmental impact; energy considerations; information systems; and design 
and engineering operations. 

1976. 706 pgs. Book No. H00097 $40.00, No member disccunt 


Symposium on Railroad Equipment Dynamics 


Railroad equipment dynamics is undergoing comprehensive study by the railroad 
industry, its supporting suppliers and the Federal Government. An understanding is 
emerging of the forcing influences involved in the dynamic environment. This Sym- 
posium contains papers involving dynamic analysis, testing of the dynamic environ- 
ment, results of freight car performance operating in the dynamic environment and a 
power system design related to improving the operational performance of a railroad 
passenger car. The papers have a permanent value for those whose investigations, 
tests and designs involve railroad vehicle vibration. 

1976. 176 pgs. Book No. G00102 $20.00, ASME members $10.00 


TRANSPORTATION TRANS 


1975 Rail Transportation Proceedings 


A collection of twenty-two technical papers covering the present state of the art of rail 
transportation. Some of the topics covered are: exhaust emissions; random rail ir- 
regularities; accidents and non-destructive inspection; retarder noise; hot spot heat- 
ing by composition shoes; aluminum center wheels for rapid transit rail service; train- 
ing center for maintenance personnel; and train derailment time-history from seis- 
mological records. 

1976. 274 pgs. Book No. GO0098 $30.00, ASME members $15.00 


TRANSPORTATION 


Other Available Proceedings 


1974 Rail Transportation Proceedings 
248 pgs. Book No. G00086 $25.00, 
ASME members $12.50 


1973 Rail Transportation Proceedings 
278 pgs. Book No. G00023 $25.00, 
ASME members $12.50 
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1972 Rail Transportation Proceedings 
296 pgs. Book No. H00053 $25.00 
ASME members $12.50 
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1971 Rail Transportation Proceedings 
228 pgs. Book No. G00080 $20.00, 
ASME members $10.00 


1970 Rail Transportation Proceedings 
220 pgs. Baok No. G00071 $19.00 
ASME members $9.50 
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